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4.  INTRODUCTION 

This  document  presents  the  research  activities  pursued  during  the  fifth  reporting  period  (June 
15,  2011  -  June  14,  2012)  of  the  project,  “ Photonic  Breast  Tomography  and  Tumor 
Aggressiveness  Assessment .”  The  major  thrust  was  on  developing  fast,  accurate,  and  noninvasive 
methods  for  detecting  and  locating  breast  tumors  in  early  growth  stages  when  those  are  more 
amenable  to  treatment.  We  continued  developing  and  testing  the  near-infrared  (NIR)  light-based 
experimental  methods  that  capitalize  on  intrinsic  differences  in  optical  absorption  and  scattering 
properties  of  normal  and  cancerous  breast  tissues.  We  have  further  initiated  the  development  of 
fluorescence  tomography  approaches  that  uses  exogenous  contrast  agents  to  enhance  contrast 
and  hold  potential  for  functional  and  molecular  imaging. 

5.  BODY 

During  this  reporting  period  we  have  built  on  our  earlier  research  on  the  project  [1-4]  and 
made  substantial  progress  in  developing  and  testing  non-invasive  near-infrared  optical  imaging 
modalities  for  early  detection  of  breast  cancer  [Specific  Aim  4\.  Our  focus  has  been  on 
developing  non-iterative  approaches  to  realize  fast,  near-real-time  image  reconstruction  of  small 
targets.  Early  detection  of  breast  tumor  with  requisite  sensitivity  and  specificity  is  a  challenging 
undertaking  and  we  continued  developing  and  testing  different  approaches.  We  have  initiated 
development  of  numerical  algorithms  for  fluorescence  tomography,  since  it  has  the  potential  to 
provide  higher  contrast  as  well  as  molecular  and  functional  information  through  the  use  of 
exogenous  fluorophores.  In  particular,  we  pursued  the  following  tasks: 

•  Comparison  of  the  efficacy  of  diffuse  optical  tomography  using  different 
Decomposition  Methods; 

•  Developing  Time  Reversal  Optical  Tomography  (TROT)  for  extended  targets;  and 

•  Developing  Fluorescence  Tomography. 

We  provide  a  brief  outline  of  activities  and  accomplishments  in  these  areas,  and  refer  to 
appended  materials  for  detailed  description  where  applicable. 

5.1.  Diffuse  Optical  Imaging  Using  Decomposition  Methods 

Diffuse  optical  imaging  (DOI)  for  detection  and  retrieval  of  location  information  of  targets  in 
a  highly  scattering  turbid  medium  may  be  treated  as  a  “blind  source  separation”  (BSS)  problem 

[5] .  BSS  is  a  general  problem  in  information  theory  that  involves  retrieval  of  “component” 
signals  from  measured  signals.  The  measured  signals  are  weighted  mixtures  of  the  component 
signals  contributed  by  the  targets  (“blind  sources”)  and  may  be  expressed  as  a  matrix  equation: 

X  =  AS,  (1) 

where  rows  of  X  represent  the  measured  mixed  signals,  rows  of  S  represent  the  component 
signals,  and  A  is  the  mixing  matrix. 

Different  matrix  decomposition  methods,  such  as,  Independent  Component  Analysis  (ICA) 

[6] ,  Principal  Component  Analysis  (PCA)  [7]  and  Non-negative  Matrix  Factorization  (NMF)  [8] 
are  used  to  solve  the  general  BSS  problem.  The  three  algorithms  have  different  assumptions, 
which  may  lead  to  different  favored  conditions.  We  introduced  these  decomposition  methods  for 
solving  the  diffuse  imaging  problem,  and  carried  out  a  systematic  study  using  simulated  and 
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experimental  data  to  compare  the  performance  of  the  three  decomposition  methods  -  ICA,  PCA, 
and  NMF  [9,  10].  We  first  reported  the  initial  results  of  this  comparative  study  in  our  fourth 
Annual  Report  [4],  The  study  involved  detection  and  retrieval  of  three-dimensional  (3-D) 
location  information  of  absorptive  and  scattering  targets  in  a  model  medium  whose  thickness 
and  optical  absorption  coefficient,  scattering  coefficient,  transport  length,  and  anisotropy  factor 
were  similar  to  the  average  values  of  those  parameters  for  a  compressed  human  breast.  The  key 
result  was  that  the  approaches  could  detect  and  extract  3-D  location  of  small  tumor-like  targets 
within  ~  2  mm  of  their  known  locations.  The  formalism  and  some  results  were  presented  in  the 
fourth  Annual  Report  [4],  and  in  the  first  part  of  this  reporting  period  we  published  a  peer- 
reviewed  journal  article  [10]  presenting  the  experimental  methods,  materials,  algorithms,  and  key 
results,  which  is  appended  (Appendix  1)  to  this  report  for  further  details. 

It  was  apparent  from  our  work  on  absorptive  and  scattering  targets  that  the  NMF  formalism, 
because  of  its  non-negativity  constraint,  may  be  particularly  suited  for  detection  of  fluorescent 
targets.  We  explored  development  of  fluorescence  tomography  (Section  5.3)  during  the  later  part 
of  the  reporting  period. 


5.2.  Time  Reversal  Optical  Tomography  (TROT)  for  Extended  Targets 

Earlier  we  reported  on  our  development  of  Time  Reversal  Optical  Tomography  (TROT)  for 
detecting  and  obtaining  3-D  location  information  about  small,  point-like  absorptive  and  scattering 
targets  [2-4,  11]  as  a  part  of  our  on-going  quest  for  fast  and  accurate  methods  for  detection  and 
localization  of  tumours  in  breast,  and  for  detection  of  margins  during  surgical  removal  of  breast 
tumours  ( Specific  Aim  4).  The  approach  combines  the  idea  of  time  reversal  (TR)  invariance  and  the 
vector  subspace  classification  method  of  Multiple  Signal  Classification  (MUSIC).  It  provided  the 
locations  of  small  absorptive  and  scattering  targets  within  a  turbid  medium  with  high  accuracy 
[11],  During  the  current  reporting  period  we  focused  on  further  developing  the  approach  to 
encompass  extended  targets  and  to  assess  their  size,  and  some  measure  of  optical  properties.  These 
assessments  are  necessary  for  target  identification,  such  as,  classifying  a  growth  in  the  breast  as 
benign  or  malignant.  We  test  the  formalism  so  developed  using  simulated  data  for  a  single 
absorptive  target,  a  single  scattering  target,  and  two  absorptive  targets  assuming  different  noise 
levels.  Here  we  present  a  brief  overview  of  the  TROT  formalism  for  small  targets  [11],  outline 
how  the  size  and  optical  strength  of  an  extended  target  may  be  estimated,  and  defer  the  details  to 
an  appended  conference  proceedings  paper  (Appendix  2). 

In  the  first  order  Born  approximation,  the  perturbation  in  the  light  intensity  distribution  due  to 
the  presence  of  the  targets  (inhomogeneities  in  optical  properties)  can  be  expressed  using  a  data 
matrix  [1 1]  of  the  form: 


IVl  1V1 

Ka \^G,(rl,Xm)raG’(Xm,r],)  =  (Xm), 


(2a) 


for  absorptive  targets,  and 

M 

dagMmym3agJ(Xm),  (2b) 

m=la={x,y,z} 

for  scattering  targets,  where  g,(r)  =  {Gs(r,  rj),  j  =  1,  . . .,  Ns)  and  gjr)  =  { Gcl(r„  r),  i=  1,  ...,  Nd  } 
are  the  Green’s  function  vectors  (GFVs)  associated  with  the  source  and  detect  planes, 
respectively;  the  superscript  T  denotes  transpose;  Gs(r,  rs)  and  Gd(rd,  r )  are  the  Green’s  functions 
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that  describe  light  propagations  in  the  background  medium  from  a  source  at  rs  to  a  target 
(inhomogeneity)  at  r  and  from  the  target  to  a  detector  at  rd,  respectively;  rm=Sjua{Xm)cSVm 
( rm  =SD{Xm)cSVm )  is  the  absorptive  (scattering)  optical  strength  of  the  mth  absorptive  (scattering) 
target  at  Xm  with  volume  SVm;  Sjua  (SD)  is  the  difference  in  the  absorption  (diffusion)  coefficients 
between  the  target  and  the  background  medium;  Ns,  Nd  and  M  are  the  numbers  of  sources, 
detectors  and  targets,  respectively.  It  is  assumed  the  number  of  targets  is  less  than  the  number  of 
sources  and  detectors,  M  <  min  (Nd,  Ns);  c  is  the  light  speed  in  the  medium.  A  time  reversal 

J.  rr t  9  rri  ^  rri 

matrix  is  constructed  as  Tsdds  =  K'K  [TDssd  =  (K  )  K  =  K  K  ]  in  frequency  domain,  and  Tsdds  = 
KrK  (TDssd  =  KKt)  in  the  continuous  wave  illumination.  TDSsd  and  Tsdds  have  a  common  set  of 
eigenvalues  j  =  1,  ...,  min(./Vs,  Nd)},  and  different  sets  of  eigenvectors  { i  =  1,  . ..,  Nd}  and 
{ vp  j  =  1,  ...,  Ns),  respectively.  The  eigenvectors  are  separated  into  signal  and  noise  subspaces 
using  an  L-curve  method  [12]  with  an  eigenvalue  threshold  e.  For  absorptive  targets,  the 
locations  are  determined  using  the  MUSIC  pseudo  spectrum  [11] 


(3a) 


associated  with  the  source  plane  or  a  similar  form  for  the  detector  plane  Pd(Xp),  or  the  product  of 
these  two, 


p(x,hp.(x,)pAx.) 


(3b) 


where  Xp  is  a  test  target  position  in  the  sample  volume.  Since  the  eigenvalues  and  eigenvectors 
of  Tsdds  and  TDssd  can  be  connected  using  singular  value  decomposition  (SVD),  i.e. 

K  ~  UYVT ,  (4) 

where  V  =  { v; }  and  U  =  { «, } ,  corresponding  to  I  =  diag{  >  e],  are  matrices  for  the  signal 

subspaces.  By  comparing  Eq.  (4)  and  Eq.  (2),  the  target  optical  property  can  be  retrieved  by 
transforming  the  eigenvalue  matrix  E  via 

r*(G‘,r1t/ZV7’((G,)rr1,  (5) 


where  T  =  diag{rm,  m  =  1,  ...,  M}\  G'  =  { gs(rm) } ,  Gd  =  { gd(rm) }  are  matrices  including  GFVs 
associated  with  the  retrieved  target  positions.  For  scattering  targets,  the  GFVs  gd  and  gs  in  Eq.(3) 
and  Eq.(5)  are  replaced  by  dagd  and  dags,  a  =  x,y,  z,  respectively. 

An  optimal  contour  (a  surface  Q  when  plotted  in  3D)  in  the  pseudo  spectrum  is  selected  to  be  the 
boundary  of  the  target(s),  via  [13] 

Q  =  argmin|^-/^(Q)|2,  (6) 

where  K  is  normalized  data  matrix  obtained  from  known  target  surface  in  simulation  (from 
experimental  measurements  for  real  targets)  and  K(Q)  is  normalized  data  matrix  calculated  from 
the  contour  of  the  pseudo  spectrum.  The  Green’s  functions  used  in  the  calculation  are  those  for 
the  intervening  medium. 

In  the  simulation  to  test  the  formalism,  the  sample  was  taken  to  be  a  40-mm  thick  uniform 
scattering  slab  with  lateral  dimension  of  80  mm  x  80  mm.  Its  absorption  and  diffusion 
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coefficients  were  assumed  to  be  pa  =  0.003  mm'1  and  D  =  1/3  mm  (transport  mean  free  path,  lt  = 
1  mm),  respectively,  which  are  similar  to  the  average  value  of  those  parameters  for  human  breast 
tissue.  The  index  of  refraction  n  of  the  medium  was  taken  to  be  1.33.  The  speed  of  light  is  2.998 
x  108  m/s,  or  299.8  mm/ns  in  vacuum,  and  225.4  mm/ns  in  the  medium.  Three  simulated  datasets 
were  generated  with  10-mm  diameter  spherical  targets  embedded.  In  the  first  dataset,  an 
absorptive  target  was  centered  at  (40,  40,  20)  mm.  In  the  second  dataset,  two  absorptive  targets 
were  located  at  (20,  40,  20)  mm  and  (60,  40,  20)  mm,  respectively.  In  the  third  dataset,  a 
scattering  target  was  centered  at  (40,  40,  20)  mm.  The  absorption  coefficient  of  all  the  absorptive 
targets  was  set  to  be  higher  than  the  background  with  A jua  =  0.001  mm'1,  while  the  diffusion 
coefficient  was  taken  to  be  the  same  as  that  of  background.  The  diffusion  coefficient  of  the 
scattering  target  was  set  to  be  lower  than  the  background  (higher  scattering  coefficient)  with  AD 
=  -0.1  mm  (lt  =  0.7  mm),  while  the  absorption  coefficient  was  taken  to  be  the  same  as  that  of  the 
background.  The  volume  of  all  targets  was  515  mm3  when  the  sample  volume  is  discretized  into 
1  mm  x  1  mm  x  1  mm  voxels  in  the  forward  model.  The  optical  strength  of  the  absorptive  targets 
was  A/uacAV  =  1 16.08  mm3/ns;  while  the  optical  strengths  of  the  scattering  target  was  ADcAV  =  - 
11608.1  mm5/ns.  The  incident  CW  beam  step  scanned  the  sample  at  41  x  41  grid  points  covering 
an  80  x  80  mm2  area,  with  a  step  size  of  2  mm.  Light  on  the  opposite  side  was  recorded  at  41  x 
41  grid  points  covering  the  same  area.  Additive  Gaussian  noise  of  different  noise  levels  was 
added  to  the  simulated  data.  The  data  matrix  K  was  then  obtained  using  Eq.  (2)  directly,  and 
analyzed  using  the  TROT  formalism.  The  results  are  shown  below.  For  simplicity,  when  the 
reconstruction  result  for  one  target  is  displayed,  a  smaller  volume  of  40  mm  x  40  mm  x  40  mm 
around  the  center  is  selected.  Due  to  the  distortion  in  the  retrieved  target  shape,  target  volume 
will  be  used  to  describe  the  target  size. 

Results  for  the  three  different  sets  of  targets  are  detailed  in  Appendix  2.  Here  we  present  the 
result  for  the  case  of  two  spherical  absorptive  targets  (diameter  10  mm)  embedded  in  the 
medium,  with  a  center-to-center  separation  of  40  mm.  Different  added  noise  levels:  0%,  5%, 
10%,  20%,  and  100%  were  considered.  The  target  locations  of  two  targets  were  accurately  found 
to  be  (20,  40,  20)  mm  and  (60,  40,  20)  mm  at  all  noise  levels.  Typical  axial  and  sagittal  images 
of  the  targets  generated  using  the  pseudo  spectra  for  a  noise  level  of  20%  are  displayed  in  Fig.  1 . 
Similar  images  were  obtained  for  other  noise  levels.  The  optical  strength  of  each  target  was 
calculated  within  3.3%  error  at  all  noise  levels.  The  size  of  targets  was  found  within  13.4%  error 
when  no  noise  was  present,  and  varied  with  noise  added.  The  retrieved  optical  strength  and  size 
of  the  targets  are  listed  in  Table  I  for  different  noise  levels.  Since  both  targets  have  the  same 
optical  property  and  size  in  the  simulated  data  and  the  retrieved  data,  only  one  set  of  values  is 
shown  in  the  table. 


axial 


_  20 

1  40 
^  60 

80 


(b)  *  (mm) 


sagittal  (x  =  20  mm)  sagittal  (x  =  60  mm) 


(c) 


(d) 


(e) 


Fig.l.  (a)  Eigenvalues  plotted  in  logarithmic  scale;  (b),  (c)  and  (d)  are  the  axial  and  sagittal  views  of  the 
targets  using  pseudo  spectrum  in  logarithmic  scale;  (e)  is  the  retrieved  target  image.  The  added  noise  level 
was  20%. 
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Further  simulations  (not  shown  here)  indicate  that  retrieved  target  size  was  more  accurate  when 
small  targets  were  involved.  Simulations  with  one  target  at  different  locations  and  two  targets 
with  different  separations  were  also  tested  for  absorptive  and  scattering  targets.  The  target 
locations  were  accurately  retrieved  in  all  of  these  cases.  The  retrieved  optical  strength  in  all  cases 
was  much  less  sensitive  to  the  target  location  and  size  than  the  retrieved  target  size.  As  expected, 
it  was  more  challenging  to  retrieve  the  optical  strength  and  size  of  scattering  targets  than 
absorptive  targets. 


Table  I.  Optical  strength  and  size  of  two  absorptive  targets  (10-mm  diameter) 


Noise  Level 
(%) 

Optical  Strength 

Size 

Retrieved 

(mm3/ns) 

Error 

(%) 

Retrieved 

(mm3) 

Error 

(%) 

0 

112.3 

3.3 

584 

13.4 

5 

112.3 

3.3 

295 

42.7 

10 

112.3 

3.3 

368 

28.5 

20 

112.3 

3.3 

487 

5.4 

100 

112.2 

3.3 

663 

28.7 

*  Known  values:  volume:  515  mm3,  optical  strength:  116.08  mm3/ns 

Further  simulations  also  showed  closer  distance  between  two  targets  made  target  size  retrieval 
more  difficult,  because  of  the  cross  talk  between  the  two  targets  showing  up  in  the  contour  of  the 
pseudo  spectrum. 

Future  work  will  involve  testing  the  efficacy  of  the  approach  using  experimental  data  for  targets 
of  different  size  and  optical  strength  with  optical  properties  similar  to  that  of  breast  tumors  at 
different  stages. 

5.3.  Fluorescence  Tomography 

Diffuse  optical  imaging  approaches  for  breast  cancer  detection  may  be  developed  capitalizing  on 
the  absorption,  scattering,  and  fluorescence  contrasts  between  the  tumor  and  the  surrounding 
tissues.  However,  there  has  been  a  recent  surge  in  interest  in  fluorescence  tomography  because  of 
its  superior  detection  sensitivity  and  specificity,  higher  signal-to-background  ratio  and  better 
spatial  resolution  than  other  diffuse  optical  imaging  (DOI)  approaches,  and  potential  to  provide 
molecular  information  on  disease-induced  changes  in  biological  tissues  [14-16],  In  addition, 
advances  towards  development  of  target-specific  exogenous  contrast  agents  [16,  17],  imaging 
instrumentation  [18,  19],  as  well  as,  analytical  methods  and  numerical  algorithms  [15,  20]  hold 
the  promise  for  noninvasive  detection  and  characterization  of  tumors.  The  use  of  target-specific 
contrast  agent  endows  fluorescence  tomography  with  potentially  higher  contrast,  and  specificity. 
The  key  attributes  of  a  clinically  useful  cancer  imaging  modality  include  early  detection, 
adequate  spatial  resolution,  speedy  image  reconstruction,  and  diagnostic  potential.  Fluorescence- 
based  approaches  are  sensitive  to  molecular  changes  and  may  provide  useful  diagnostic 
information,  but  diffuse  nature  of  light  propagation  in  biological  tissues  impedes  spatial 
resolution.  For  these  reasons,  we  initiated  the  development  of  fluorescence  tomography  towards 
the  end  of  the  current  reporting  period. 

We  plan  to  develop  fluorescence  tomography  approaches  based  on  the  TROT  and  NMF,  which 
we  previously  explored  and  tested  for  absorptive  and  scattering  targets  [11,  9].  While  other 
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fluorescence  tomography  approaches  are  computation  intensive  since  those  try  to  iteratively 
compute  fluorescence  properties  (such  as,  fluorophore  concentration)  of  every  voxel  of  the 
sample  volume,  both  TROT  and  NMF  are  non-iterative  and  are  expected  to  provide  faster 
localization.  During  this  reporting  period  we  started  developing  the  algorithm  for  fluorescence 
TROT  and  used  it  with  experimental  data  for  detecting  and  locating  a  single  fluorescent  target  in 
a  breast  simulating  turbid  medium. 

The  experimental  arrangement  uses  a  multi-source  illumination  and  multi-detector  signal 
acquisition  scheme  to  acquire  multiple  angular  views  of  the  sample.  The  input  surface  (source 
plane)  of  the  slab  sample  is  scanned  across  an  external  point  light  source  (laser  beam).  The 
fluorescent  target(s)  embedded  in  the  highly  scattering  medium  are  excited  by  diffusely 
propagating  light  of  wavelength  Xx.  A  fraction  of  the  forward  propagating  fluorescence  signal 
emitted  by  the  targets  at  wavelength  Xm,  is  detected  on  the  other  side  of  the  sample  by  a  two- 
dimensional  detector  array.  For  small  fluorescent  targets  at  r;  with  volume  V),  the  fluorescence 
signal  from  the  target  at  r  illuminated  by  a  point  source  of  unit  power  at  rs  is  given  by[20] 

K  =  (r,  ’  0J)fj  (oj)S,T  (rj  ,#>),  (7) 

j 

where  gs  (r,  to)  =  {Gx  ( r ,  r„  co)}1  and  gd  (r,  to)  =  {Gm  (rd,  r,  co)}T,  (the  superscript  T  denotes 
transpose);  Gx  (r,  rs,  co)  is  a  Green’s  function  that  describes  the  propagation  of  excitation  light  at 
excitation  wavelength  Xx  from  the  source  at  rs  to  the  target  at  r;  Gm  (rd,  r,  co)  is  a  Green’s  function 
that  describes  the  propagation  of  the  fluorescent  light  at  emission  wavelength  Xm  from  the  target 
at  r  to  the  detector  at  rd;  co  is  the  modulation  angular  frequency  of  the  light;  f/co)  is  the 
fluorescence  strength  of  the  jth  target,  given  by 

fj  M  =  yirj  )cmVj  /  [1  -  icoT(rj )] , 


y  is  the  fluorescence  yield,  cm  is  the  speed  of  light  in  the  medium,  and  r  is  the  fluorescence 
lifetime.  K  describes  that  the  diffuse  propagation  of  excitation  light  of  wavelength,  Xx  from  the 
sources  through  the  medium  to  illuminate  the  targets,  and  then  the  propagation  of  emitted 
fluorescence  of  wavelength,  Xm  from  the  targets  to  the  detectors.  In  this  study,  continuous  wave 
(CW)  illumination  is  used,  i.e.  co  =  0.  A  time  reversal  matrix  Tsdds  =  KTK  ( Tdssd  =  KKT)  is 
constructed.  A  set  of  eigenvectors  {uk,  k  =  1,  ...,  Nd)  and  [v/,  /  =  1,  ...,  Ns}  are  calculated  for 
Tdssd  and  Tsdds ,  respectively,  with  common  eigenvalues  {jUj,  j  =  1,  ...,  min(Ai?  Nd) } ,  where  Ns 
and  Nd  are  numbers  of  sources  and  detectors,  respectively[ll].  The  eigenvalues 

jUj  =  |  /)  | 2  jgd  (r7 ,  <w)|2  ||gv  (/-y ,  are  proportional  to  squared  fluorescence  strengths  of  the  targets  if 

the  targets  are  well  resolved;  otherwise,  they  are  linear  combinations  of  individual  fluorescence 
strengths  of  different  targets  [11]. 


By  using  an  L-curve  method  with  an  eigenvalue  threshold  s  [12],  eigenvectors  are  separated  into 
signal  and  noise  subspaces.  A  pseudo  spectrum  associated  with  detector  plane  is  calculated  using 
multiple  signal  classification  (MUSIC)  for  a  test  target  position  Xp  in  the  sample  volume  [1 1] 


pAxp^)=1  I 


jUj<£ 


sMp’V) 


8d{Xp ’®) 


(9) 


The  locations  of  targets  are  retrieved  to  be  the  poles  of  the  pseudo  spectrum.  A  similar  pseudo 
spectrum  for  the  source  plane, 
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may  also  be  used  to  retrieve  the  target  position. 

The  sample  used  in  the  experiment  to  demonstrate  the  approach  was  a  250  mm  x  250  mm  x  60 
mm  rectangular  transparent  plastic  cell  filled  with  Intralipid-20%  (Baxter)  suspension  in  distilled 
water  with  a  fluorescent  targets  embedded  inside.  The  thickness  of  the  sample  was  comparable  to 
that  of  a  typical  compressed  breast.  The  concentration  of  Intralipid-20%  was  adjusted[21]  to 
provide  a  transport  mean  free  path  lt  of  ~lmm,  which  happens  to  be  about  the  same  for  both 
excitation  and  emission  wavelengths,  and  was  similar  to  the  average  value  of  lt  for  human  breast 
tissue  at  these  wavelengths.  The  fluorescent  target  was  a  4.2  mm  inner-diameter  10-mm  long 
cylindrical  glass  tube  filled  with  solution  of  ICG  dye  (Sigma- Aldrich)  water.  The  dye  solution  in 
the  targets  was  prepared  by  dissolving  ICG  at  a  concentration  of  1  pM  in  the  Intralipid-20% 
suspension  of  same  concentration  as  the  background  medium  to  ensure  that  the  target  had  the 
same  scattering  coefficient  as  the  background  medium,  but  a  higher  absorption  coefficient  of 
0.027  mm1.  The  water  solution  of  ICG  absorbs  light  over  the  600  -  900  nm  range  with  peak  at 
780  and  fluoresces  in  the  790  -  966  nm  range  with  peak  at  around  825  nm  in  the  NIR  enabling 
deeper  penetration  of  light  in  tissues.  It  has  been  approved  by  Food  and  Drug  Administration 
(FDA)  for  biomedical  applications. 

The  experimental  arrangement  is  shown  schematically  in  Fig.  2  ( Specific  Aim  4,  Task  14).  The 
target  was  embedded  in  the  mid-plane  ( z  =  30  mm)  of  the  sample  cell.  A  multi-source  sample 
excitation  and  multi-detector  fluorescence  signal  acquisition  scheme  was  used  to  acquire 
multiple  angular  views  of  the  sample. 


(a) 


Laser 


Sample 


x 

Source 

Plane 


Raw  Image  Cropped  Image  Binned  Image 


Fig. 2.  (a)  A  schematic  diagram  of  the  experimental  arrangement  for  imaging  objects  embedded  in  a  turbid 
medium.  [Key:  NB  =  narrow  band  pass  filter,  TS  =  translational  stage,  CCD  =  charge  coupled  device,  PC  = 
personal  computer]  Inset  (below)  shows  the  2-D  array  in  the  input  plane  that  was  scanned  across  the 
incident  laser  beam  and  a  typical  raw  image  is  shown  in  the  PC  monitor,  (b)  The  absorption  and 
fluorescence  spectra  of  ICG  in  water,  (c)  A  typical  raw  image  is  cropped  and  binned. 
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The  entrance  face  of  the  slab  sample  (source  plane)  was  illuminated  by  a  100-mW  790-nm  diode 
laser  beam.  The  multi-source  illumination  scheme  was  realized  by  scanning  the  sample  across  the 
laser  beam  in  a  two-dimensional  x-y  array  of  9  x  9  grid  points  with  a  5-mm  step  size  using  a 
computer-controlled  translation  stage.  A  cooled  CCD  camera  equipped  with  a  60-mm  focal-length 
camera  lens  collected  and  sensed  a  fraction  of  the  forward  propagating  fluorescence  signal  from 
the  opposite  face  of  the  sample  (detection  plane )  through  a  narrow-band  (full-width-at-  half¬ 
maximum  (FWHM)  bandwidth  of  10  nm)  interference  filter  with  50%  peak  transmission  at  830 
nm.  The  CCD  camera  had  1024  x  1024  pixels  with  a  pixel  size  of  24  pm,  and  was  considered  as 
multi-detector,  since  each  illuminated  pixel  could  be  considered  as  an  individual  detector.  The  NB 
filter  was  chosen  to  select  a  substantial  fraction  of  the  ICG  fluorescence  around  the  peak  emission 
wavelength  and  block  the  transmitted  790  nm  pump  light.  The  scanning  and  data  acquisition 
processes  were  controlled  by  a  personal  computer  (PC).  Raw  images  were  recorded  by  the  PC  for 
each  scan  position,  and  stored  for  subsequent  analysis.  A  typical  image,  which  is  a  2 -D  intensity 
distribution,  is  shown  in  the  left  frame  of  Fig.  2(c). 

At  every  scan  position,  an  image  was  also  acquired  with  the  NB  filter  removed.  In  this  case,  the 
recorded  images  were  essentially  transillumination  images  since  the  fluorescence  signal  was 
negligible  compared  to  the  much  more  intense  transmission  signal  (ratio  of  transmission  signal  to 
fluorescence  signal  -1500).  Thus  correlated  imaging  and  retrieval  of  target  location  using  both 
fluorescence  and  transmission  measurements  were  enabled  in  this  experiment.  The  transmission 
images  were  analyzed  to  estimate  the  average  value  of  the  optical  property  of  the  medium  k  = 
(3paMs’)V2,  where  pa  and  p, ’  are  the  absorption  and  reduced  scattering  coefficients  at  790  nm, 
respectively.  In  fact,  the  values  of  these  optical  parameters  of  Intralipid-20%  suspension  in  water 
happened  to  be  very  close  for  excitation  and  fluorescence  wavelengths. 

From  each  fluorescence  image,  a  region  of  interest  was  cropped  out  and  every  5x5  pixels  in  the 
cropped  image  were  binned  to  one  pixel  to  enhance  the  signal-to-noise  ratio.  The  response  data 
matrix  was  constructed  using  the  transmitted  fluorescence  light  intensity  distribution  in  the 
processed  images.  The  TR  matrix  was  generated  by  multiplying  the  response  matrix  by  its 
transpose  for  our  CW  probing  scheme.  The  eigenvalue  equation  of  TR  matrix  was  solved.  Then  the 
signal  and  noise  subspaces  were  separated.  MUSIC  pseudo  spectrum  was  calculated  and  target 
locations  were  determined  using  the  poles  in  the  pseudo  spectrum. 

For  comparison,  the  transmission  data  were  also  analyzed  using  TROT  as  detailed  in  Ref.  11.  In 
this  case,  the  target  was  absorptive,  and  the  contrast  was  mainly  due  to  higher  absorption  of  the 
excitation  beam  by  the  target.  It  should  be  noted  that  absorption  measurement  involves  changes  in 
the  intensity  of  the  excitation  beam,  and  consequently  the  TROT  analysis  used  the  difference 
images  between  the  raw  transmission  images  and  a  reference  image  for  the  background  medium. 
The  reference  image  may  be  estimated  as  the  average  of  the  images  acquired  at  all  scan  positions. 

The  TR  matrix  was  constructed  using  the  fluorescence  data,  and  an  eigenvalue  equation  was  then 
solved.  The  eigenvector  with  the  dominant  eigenvalue  was  used  to  calculate  the  pseudo  spectrum 
for  all  voxels  in  the  3-D  space  of  the  sample.  The  voxel  size  was  0.77  mm  x  0.77  mm  x  1  mm. 
Three-dimensional  tomographic  pseudo  images  were  generated  using  the  pseudo  spectrum.  The 
single  target  was  detected,  and  the  position  of  the  target  was  determined  using  the  peak  in  the 
pseudo  spectrum  and  listed  in  Table  II,  with  comparison  to  the  actual  position.  The  image  at  the 
retrieved  z-coordinate  of  the  target  position  (z  =  30.5  mm)  plotted  using  the  pseudo  spectrum  is 
shown  in  Fig.  3(a). 
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Fig.  3  TROT-reconstructed  image  at  z  =  30.5  mm  using  fluorescence  data  is  shown  in  (a),  and  at  z 
=  29.5.mm  using  transmission  data  shown  in  (b). 


Table  II.  Known  and  retrieved  target  positions. 


Known 

Retrieved 

Fluorescence 

Transmission 

Position  [x,  y,  z  (mm)] 

25.1,25.7,  30.0 

24.1,25.6,  30.5 

21.8,  25.6,  29.5 

Error  [Ax,  Ay,  Az  (mm)] 

- 

1.0,  0.1,  0.5 

3.3,  0.1,  0.5 

FWHM  [Sc,  Sy  (mm)] 

- 

3.8,  3.8 

6.9,  6.9 

The  transmission  data  was  then  analyzed  for  comparison.  The  target  was  also  detected  and  its 
retrieved  location  is  listed  in  Table  I  for  comparison.  The  image  of  the  target  at  z  =  29.5  mm  is 
shown  in  Fig.  3(b). 

As  shown  in  Table  II,  the  location  of  target  retrieved  from  the  fluorescence  data  is  in  excellent 
agreement  with  the  known  position,  and  is  consistent  with  that  retrieved  using  transmission  data. 
The  pole  of  the  pseudo  image  using  fluorescence  data  is  sharper  than  that  obtained  using 
transmission  data.  The  FWHM  of  the  pole  in  both  x  and  y  directions  in  the  fluorescence-TROT 
image  is  3.8  mm  and  6.9  mm  in  the  transmission-TROT  image. 

The  results  show  that  fluorescence  TROT  approach  could  detect  and  retrieve  the  locations  of  a 
small  fluorescent  target  embedded  in  a  breast-simulating  turbid  medium  with  the  thickness 
comparable  to  that  of  a  realistic  compressed  breast.  The  location  of  targets  was  retrieved  with  an 
accuracy  of  ~1  mm  in  all  three  dimensions.  We  plan  to  pursue  this  further  for  multiple  targets. 
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5.4.  Research  Proposal  Development 

One  of  the  proposed  tasks  ( Specific  Aim  0,  Task #  7)  involves  developing  at  least  one  research 
proposal  and  submitting  it  to  NIH  or  USAMRMC  for  funding.  We  submitted  a  pre-proposal 
entitled,  “ Multimodal  nanocomposites  for  detection  and  prevention  of  breast  cancer  metastases ” 
to  the  Idea  Award  (Collaborative  Option)  category  of  the  2012  Breast  Cancer  Research  Program 
of  CDMRP.  S.  K.  Gayen  (PI  of  this  proposal)  is  the  Initiating  PI  of  the  above-mentioned 
proposal.  The  Collaborative  PI  is  Valeria  Balogh-Nair  of  the  CCNY  Chemistry  Department.  A 
brief  overview  of  the  pre-proposal  follows. 

Research  Idea 

The  objective  of  the  proposed  research  is  to  develop  multi-functional,  tumor-targeting 
nanocomposites  and  near-infrared  (NIR)  optical  imaging  approaches  for  early  detection  of  breast 
cancer,  and  prevention  of  micro-metastases  that  are  responsible  for  majority  of  breast  cancer 
mortality.  The  nanocomposite  synthesis  will  use  a  multivalent  dendrimer  (a  class  of  organic 
macromolecules)  platform  to  incorporate  fluorescent  moieties  ( e.g .,  semiconductor  quantum 
dots,  or  gold/silver  nanoparticles)  as  contrast  agents  for  imaging ,  and  chemokine  mimics  for 
selective  targeting  of  cancer  cells  and  prevention  of  metastases.  The  high  affinity  of  chemokine 
mimic’s  ligands  for  chemokine  CXCR4  receptors  will  ensure  selective  delivery  of  the 
nanocomposite  to  the  target.  A  multi-source  NIR  probing  and  multi-detector  signal  acquisition 
arrangement  along  with  a  time-reversal  image  reconstruction  algorithm  will  be  used  for  fast 
detection  of  tumor  and  determination  of  its  location  in  three  dimensions. 

Impact  and  Innovation 

The  major  impact  of  the  proposed  research  is  that  it  has  the  potential  to  provide  a  modality  for 
early  detection  of  breast  cancer,  and  for  prevention  of  metastases,  the  major  cause  of  mortality.  A 
broader  potential  impact  is  that  the  approach  could  be  extended  to  other  cancers  using 
chemokine  mimics  directed  to  other  chemokine  receptors. 

The  project  is  highly  innovative  in  many  ways.  First,  the  synthesis  process  brings  together  the 
concept  of  multivalency,  the  idea  of  chemokine  mimics  for  prevention  of  cancer  cell  migration, 
and  use  of  dendrimers  as  stabilizers  and  nanoreactors  for  contrast  agent  synthesis.  A  combination 
of  these  novel  ideas  is  proposed,  to  the  best  of  our  knowledge  for  the  first  time,  as  a  multi-prong 
approach  to  fight  the  menace  of  breast  cancer.  Second,  the  affinity  of  the  chemokine  mimics  to 
seven  helix  chemokine  CXCR4  receptors  obviates  the  need  for  a  targeting  vector.  The 
chemokine  mimics  will  play  the  dual  role  of:  (1)  “homing  devices”  for  selective  delivery  of  the 
nanocomposite  to  the  tumor  sites;  and  (2)  “prevention  agents”  interacting  with  the  chemokine 
receptor  sites  to  inhibit  metastases.  Third,  the  efficacy  of  noninvasive  NIR  imaging  approach  for 
early  detection  and  potential  diagnosis  will  be  significantly  enhanced  through  design  of  efficient 
contrast  agent.  Fourth,  the  idea  of  time-reversal  optical  tomography  (TROT)  is  a  new  paradigm 
in  diffusive  optical  tomographic  (DOT)  imaging.  While  currently  pursued  DOT  approaches  are 
iterative  and  computation  time  intensive,  TROT  is  non-iterative  and  will  be  faster,  which  is  a 
necessary  condition  for  real-time  imaging.  TROT  is  designed  for  detecting  and  locating  small 
targets  and  will  be  suited  for  early  detection  when  tumors  are  small.  Finally,  the  dendrimer 
platform  with  multivalent  surface  will  enable  incorporation  of  moieties  that  enhance  other 
imaging  modalities,  such  as,  magnetic  resonance  imaging,  enabling  development  of  sought-after 
dual  or  multimodal  imaging  modules. 
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6.  KEY  ACCOMPLISHMENTS 

Key  accomplishments  include: 

(a)  Development  of  a  diffuse  optical  imaging  approach  based  on  decomposition  methods 
(ICA,  PCA  and  NMF),  comparison  of  the  performance  of  the  three  algorithms  using 
experimental  data  on  small  tumor  simulating  targets  embedded  in  breast  simulating 
model  media  sample,  and  publication  of  the  results  in  a  peer-reviewed  journal  article; 

(b)  Extension  of  Time  Reversal  Optical  Tomography  (TROT)  (which  is  fast  and  minimally 
iterative)  to  extended  targets  and  testing  of  its  efficacy  to  locate  the  target  and  assess  its 
size  and  optical  strength  using  simulated  data; 

(c)  Initiating  the  development  of  the  TROT  formalism  for  fluorescent  targets,  and  testing  it 
for  a  single  target  in  a  breast  phantom  composed  of  a  model  medium  that  simulates  the 
size  and  average  optical  properties  of  human  breast; 

(d)  Presentation  of  research  results  at  the  SPIE's  International  Symposium  on  Biomedical 
Optics,  BiOS  '12/  Photonics  West,  21-26  January,  2012,  San  Francisco,  California. 

(e)  Recognition  of  the  TROT  formalism  through  its  highlighting  in  “ Optics  in  2012 ”  a 
special  of  the  Optics  and  Photonics  News  (published  by  the  Optical  Society  of  America, 
December  2012).  The  OPN  noted,  "Bob  Guenther,  the  guest  editor  of  the  issue,  and  a 
team  of  seven  volunteer  editors  have  combed  through  the  work  of  scientists  from  around 
the  globe  to  identify  30  summaries  that  highlight  the  most  exciting  optics  research  to 
emerge  over  the  preceding  12  months.”  The  journal  article  presenting  our  work  on  TROT 
[“Time  reversal  optical  tomography:  locating  targets  in  a  highly  scattering  turbid 
medium,”  Opt.  Express  19,  21956-21976  (2011)]  was  summarized  for  this  special  issue. 

(f)  Development  and  submission  of  an  independent  research  pre-proposal  to  BCRP  2012,  an 
indication  of  the  CCNY  team’s  progress  towards  developing  a  sustainable  breast  cancer 
research  program  at  CCNY. 

7.  REPORTABLE  OUTCOMES 

Publications 

(1)  Binlin  Wu,  M.  Alrubaiee,  W.  Cai,  M.  Xu,  and  S.  K.  Gayen,  “Diffuse  optical  imaging 
using  decomposition  methods,”  International  J.  Opt.  2012,  Article  ID  185435,  12  pages 
(2012);  doi:10. 1155/2012/185435.  {Attached  as  Appendix  1) 

(2)  Binlin  Wu,  W.  Cai,  M.  Xu,  and  S.  K.  Gayen,  “Time-reversal  optical  tomography: 
detecting  and  locating  extended  targets  in  a  turbid  medium,”  Proc.  SPIE  Vol.  8216 
82160K-1  -  K-5  (2012),  doi:  10.1117/12.906799  {Attached  as  Appendix  2) 

(3)  Binlin  Wu,  Wei  Cai,  Swapan  K.  Gayen,  “Time  reversal  optical  tomography,”  Optics  and 
Photonics  News  (OPN)  special  issue  “Optics  in  2012”,  December  2012,  p.  41.  Although 
the  summary  was  published  in  the  December  2012  issue  of  OPN,  the  source  paper  was 
published  in  2011.  So,  we  include  it  in  this  late  report.  {Attached  as  Appendix  3) 
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Conference  Presentations 

Binlin  Wu,  W.  Cai,  M.  Alrubaiee,  M.  Xu  and  S.  K.  Gayen,  “Time-reversal  optical  tomography: 
detecting  and  locating  extended  targets  in  a  turbid  medium,”  Paper  8216-19  presented  at  the 
SPIE's  International  Symposium  on  Biomedical  Optics,  BiOS  '12/  Photonics  West,  21-26 
January,  2012,  San  Francisco,  California. 

Research  Proposal 

S.  K.  Gayen  (Initiating  PI),  Y.  Balogh-Nair  (Collaborating  PI),  “ Multimodal  nanocomposites 
for  detection  and  prevention  of  breast  cancer  metastases submitted  to  the  Idea  Award 
(Collaborative  Option)  category  of  the  2012  Breast  Cancer  Research  Program  of  CDMRP. 

8.  CONCLUSION 

The  work  carried  out  during  this  reporting  period:  (a)  shows  the  potential  for  noninvasive 
detection  and  three-dimensional  localization  of  a  tumor  within  a  breast  with  significant  accuracy 
based  on  the  differences  in  the  light  scattering  and  absorption  characteristics  of  the  tumor  and 
normal  breast  tissue;  (b)  presents  the  formalism  and  initial  results  of  a  fluorescence  tomography 
approach  that  may  provide  higher  contrast  and  better  target  identification  based  on  molecular 
signatures  and  target  selectivity  of  exogenous  contrast  agents. 

“So  What  Section,, 

•  The  National  Cancer  Institute  (NCI)  has  identified  the  development  of  imaging 
methodologies  as  an  extraordinary  opportunity  for  advancement  in  cancer  research.  Since  the 
background  of  the  CCNY  team  is  in  physical  sciences  and  engineering,  the  training  they 
received  has  provided  them  with  necessary  laboratory  background  in  the  biology  of  cancer 
research,  and  helping  develop  a  knowledgeable  multidisciplinary  research  force  in  the  fight 
against  breast  cancer. 

•  A  recent  study  involving  35,319  patients  underscores  the  influence  of  primary  tumor  location 
on  breast  cancer  prognosis  [22],  and  makes  it  imperative  that  breast  cancer  detection 
modalities  to  obtain  three-dimensional  (3-D)  location  of  the  tumor  relative  to  the  axilla  be 
developed.  The  optical  imaging  techniques  (TROT,  fluorescence  TROT  and  other 
decomposition  methods)  being  developed  are  an  important  steps  in  obtaining  3-D  location  of 
a  tumor  within  the  breast.  The  methods  are  minimally  iterative,  fast,  and  designed  for 
locating  small  targets,  which  make  those  suitable  for  detecting  tumors  in  early  stages  of 
development  when  those  are  more  amenable  to  treatment. 
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Diffuse  optical  imaging  (DOI)  for  detecting  and  locating  targets  in  a  highly  scattering  turbid  medium  is  treated  as  a  blind  source 
separation  (BSS)  problem.  Three  matrix  decomposition  methods,  independent  component  analysis  (ICA),  principal  component 
analysis  (PCA),  and  nonnegative  matrix  factorization  (NMF)  were  used  to  study  the  DOI  problem.  The  efficacy  of  resulting 
approaches  was  evaluated  and  compared  using  simulated  and  experimental  data.  Samples  used  in  the  experiments  included 
Intralipid- 10%  or  Intralipid-20%  suspension  in  water  as  the  medium  with  absorptive  or  scattering  targets  embedded. 


1.  Introduction 

Diffuse  optical  imaging  (DOI)  for  detection  and  retrieval  of 
location  information  of  targets  in  a  highly  scattering  turbid 
medium  may  be  treated  as  a  blind  source  separation  (BSS) 
problem  [1,  2].  Various  matrix  decomposition  methods, 
such  as,  independent  component  analysis  (ICA)  [3],  princi¬ 
pal  component  analysis  (PCA)  [4],  and  nonnegative  matrix 
factorization  (NMF)  [5,  6]  have  been  developed  for  solving 
the  BSS  problem  and  retrieving  desired  information. 

Xu  et  al.  adapted  ICA  of  information  theory  to  develop 
optical  tomography  using  independent  component  analysis 
(OPTICA)  and  demonstrated  its  application  for  diffuse 
imaging  of  absorptive,  scattering,  and  fluorescent  targets 
[7-11].  ICA  assumes  the  signals  from  different  targets  to 
be  independent  of  each  other  and  optimizes  a  relevant 
measure  of  independence  to  obtain  the  ICs  associated  with 
different  targets.  The  position  coordinates  of  targets  in  three 
dimensions  are  determined  from  the  individual  components 
separately. 

PCA  assumes  that  the  PCs  contributing  to  the  signal  are 
uncorrelated  and  explain  the  most  variance  in  the  signal. 
PCA  has  been  widely  used  in  various  applications,  such  as 
spectroscopy  [12],  face  recognition  [13],  and  neuroimaging 
[14].  NMF  seeks  to  factorize  a  matrix  into  two  nonnegative 
matrices  (component  signals  and  weights)  and  requires  the 


contributions  to  signal  and  the  weights  of  the  components 
to  be  non-negative.  It  does  not  imply  any  relationship 
between  the  components.  NMF  has  also  been  widely  used  in 
biological  analysis  [15]  and  spectral  analysis  [16]. 

The  objective  of  this  study  is  to  test  and  compare 
the  efficacy  of  these  algorithms  when  used  to  solve  the 
DOI  problem.  Results  are  presented  and  compared  using 
simulative  data  and  experimental  data  using  absorptive  and 
scattering  targets  embedded  in  model  scattering  media. 
Our  interest  in  solving  the  DOI  problem  derives  from  the 
need  for  a  noninvasive  modality  for  detecting,  locating,  and 
diagnosing  breast  tumors  in  early  stages  of  growth. 

The  remainder  of  the  paper  is  organized  as  follows.  In 
Section  2,  the  formalisms  of  the  three  methods  are  intro¬ 
duced.  Section  3  evaluates  the  resulting  imaging  approaches 
using  simulated  data.  The  approaches  are  further  examined 
in  Section  4  for  experimental  data  acquired  using  absorptive 
and  scattering  targets  embedded  in  model  scattering  media. 
Section  5  summarizes  and  discusses  the  results. 

2.  Formalism 

2.1.  Blind  Source  Separation  Problem.  Blind  source  sepa¬ 
ration  (BSS),  also  known  as  blind  signal  separation,  is  a 
general  problem  in  information  theory  that  seeks  to  separate 
different  individual  signals  from  the  measured  signals,  which 
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are  weighted  mixtures  of  those  individual  signals.  Assuming 
M  individual  signals,  Sj(t),  j  =  1 .  .M,  are  linearly  mixed 
instantaneously,  the  BSS  problem  is  modeled  as  follows.  The 
dimension  of  Sj(t)  is  Ns,  the  number  of  sampling  times. 
In  this  study,  t  will  be  replaced  by  spatial  positions  of  the 
excitation  light  sources.  A  total  of  Nd  detectors  sense  Nd 
different  mixtures  of  Sj(t).  The  mixture  measured  by  the 
zth  detector  can  be  presented  as  xi(t)  =  Xjli  aySj(t),  or 
X  =  AS,  in  a  matrix  notation,  where  A  E  RNdxM  is  a 
mixing  or  weighting  matrix,  S  E  RMxN*y  X  E  RNdxN* ,  and 
M  <  min (Ns,Nd).  The  objective  of  BSS  is  to  retrieve  the 
signals  Sj(t)  and  their  weights,  ay.  ICA,  PCA,  and  NMF  are 
statistical  analysis  methods  used  to  solve  the  BSS  problem. 

2.2.  Diffuse  Optical  Imaging  Problem.  In  DOI,  one  measures 
the  signal  at  the  sample  boundary,  which  includes  a  weighted 
mixture  of  contributions  from  embedded  targets.  One  uses 
the  diffusion  approximation  [17-19]  of  the  radiative  transfer 
equation  [20,  21]  as  the  forward  model  to  describe  light 
propagation  in  a  highly  scattering  turbid  medium.  The 
perturbation  in  the  light  intensity  distribution  measured 
on  the  boundary  of  the  sample  due  to  the  presence  of  the 
targets  (which  are  localized  inhomogeneities  in  the  optical 
properties  within  the  sample  volume)  may  be  written,  in  the 
first-order  Born  approximation,  as  [22,  23] 

A0(rd,rs)  =  -  G(rd,r)8[ia(r)cG(r,rs)d3r 

J  (1) 
<5D(r)cVrG(rd,r)  ■  VrG(r,rs)d3r, 

where  r5,  r j,  and  r  are  the  positions  of  a  source  of  unit  power, 
detector  and  target,  respectively;  G(r,  r5)  and  G(rj,  r)  are  the 
Green  s  functions  that  describe  light  propagation  from  the 
source  to  the  target  and  from  the  target  to  the  detector, 
respectively;  8pa  and  8D  are  the  differences  in  absorption 
coefficient  and  diffusion  coefficient  between  the  targets  and 
the  background  medium,  respectively;  and  c  is  the  light  speed 
in  the  medium. 

A  multisource  illumination  and  multidetector  signal 
acquisition  scheme  is  used  to  acquire  light  transmitted 
through  a  scattering  medium.  For  small  absorptive  targets, 
a  perturbation  data  matrix  is  constructed  using  -A 0  for 
all  sources.  The  elements  of  the  data  matrix  pertaining  to 
absorptive  targets  represented  by  the  first  term  in  (1)  may 
be  written  in  a  discrete  form  as 

M 

Xij  =  X  Gd(ri,rm)rmGs(rm,r;) 

m=  1  V2) 

(i=  l,2,...,Nd;j  =  1,2,...  ,Ns), 

where  riy  r;-,  and  rm  are  the  locations  of  the  zth  detector, 
jth  source  and  mth  target,  respectively;  Ns,  Nd ,  and  M  are 
the  numbers  of  sources,  detectors,  and  targets,  respectively; 
rm  =  8pa(rm)c8Vm  is  the  optical  absorption  strength  of  the 
mth  target  of  volume  8Vm ;  G5(rm,r;  )  and  Gd(r/,rm)  are  the 
Greens  functions  that  describe  light  propagation  from  jth 
source  to  mth  target  and  from  mth  target  to  zth  detector, 


respectively.  The  number  of  targets  is  assumed  to  be  less  than 
that  of  sources  and  detectors,  M  <  min (Nd,Ns). 

The  mth  target  may  be  considered  to  be  a  virtual  source 
of  strength  rmG5(rm,ry)  excited  by  the  real  light  source 
located  at  r  j.  The  data  matrix  X  =  {Xy}  may  be  considered 
to  be  a  set  of  combinations  of  light  signals  from  all  virtual 
sources  mixed  by  a  mixing  matrix  {Grf(r*,rm)}.  Therefore, 
this  problem  can  be  treated  as  a  BSS  problem. 

As  the  second  term  in  ( 1 )  suggests,  each  scattering  target 
is  represented  by  three  colocated  virtual  sources  of  strength: 
rmdpG5(rm,r;),  where  dp  =  d/dp ,  (p  =  x,y,z),  and  rm  = 
8D(rm)c8Vm,  is  the  optical  scattering  strength  of  the  mth 
target  [8].  The  mixing  matrices  become  {dpGd(ri,  rm)}  for 
the  three  virtual  sources  generated  by  the  mth  target.  The 
elements  of  the  data  matrix  for  scattering  targets  may  be 
written  as 

M 

Xij  =  X  X  apGJ(ri,rm)rmapGs(rm,r;).  (3) 

m= 1  P={x,y,z} 

Since  one  absorptive  target  is  represented  by  one  cen- 
tro symmetric  virtual  source,  while  three  virtual  sources  (one 
centrosymmetric  and  two  dumb-bell  shaped)  represent  one 
scattering  target  [7,  8],  the  number  and  patterns  of  virtual 
sources  may  be  used,  in  favorable  situations,  to  identify  the 
target  as  absorptive  or  scattering  in  nature.  In  this  paper,  only 
small  targets  are  considered  since  all  three  algorithms  are 
suited  for  small  targets,  and  early  detection,  when  the  tumors 
are  more  amenable  to  treatment,  is  of  practical  interest. 

2.3.  DOI  as  a  BSS  Problem.  The  data  matrix  for  the  DOI 
problem  may  be  written  as 

M 

X  =  AS  =  XAimSmj,  (4) 

m=  1 

where  A  E  £^xM,  S  e  RMxNs ,  and  X  E  RN^N>m  For 
absorptive  targets, 


Aim  —  fim G  (r;,rm),  Smj  —  amG  ^rm,r (3a) 
while  for  scattering  targets, 

Aim  =  PmdpG  (r/,rm),  Smj  =  ocmdpG  (3b) 

{Smj}  ( j  =  1,2, ...  ,Ns  )  and  {Aim}  (f  =  1,2,.  ..Nd)  are 
two-dimensional  intensity  distributions  on  the  source  and 
detector  planes,  respectively.  Source  and  detector  planes  are 
the  boundaries  of  the  sample  through  which  light  enters  and 
exits  the  sample  volume,  respectively.  The  scaling  factors 
and  am  are  related  to  the  target  optical  strength,  rm  =  ocml3m- 
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The  location  of  the  target  and  the  scaling  factors  can  be 
retrieved  using  a  least  squares  fitting  via 


for  absorptive  and  scattering  targets,  respectively.  However, 
when  a  scattering  target  is  embedded  deep  in  a  turbid 
medium,  only  the  rmdzG5(rm,r;)  virtual  source  remains 
significant.  So,  only  p  =  z  may  be  used  for  fitting  in  (6b) 
[8]. 

2.3.1.  ICA.  OPTICA  assumes  that  the  virtual  sources  are 
independent  of  each  other  [8].  So,  they  can  be  retrieved 
through  an  iterative  process  which  seeks  to  maximize  the 
independence  among  the  components.  In  practice,  the  inde¬ 
pendent  components  are  found  by  maximizing  some  mea¬ 
sure  of  non-Gaussianity,  such  as  kurtosis  (the  fourth- order 
cumulant),  of  the  unmixed  components.  A  Matlab  program 
for  ICA  was  adopted  from  http://sccn.ucsd.edu/eeglab/. 
The  location  of  the  target  can  be  retrieved  by  fitting  the 
independent  component  intensity  distributions  (ICIDs)  to 
Green  s  functions  or  derivatives  of  Green  s  functions  using 
(6a)  and  (6b). 

2.3.2.  PC  A.  PCA  assumes  that  the  virtual  sources  are 
uncorrelated  so  that  the  correlation  (covariance)  between 
them  is  ideally  zero  and  minimal  in  practice.  The  covariance 
matrix  of  S,  cov(S)  should  be  diagonal.  The  general  process 
of  PCA  is  as  follows.  The  data  matrix  X  =  AS  +  W,  where 
Af  is  random  noise  added  to  the  data,  A  and  S  the  same  as 
defined  in  (4).  When  S  is  mean  centered,  elements  of  the 
mean- centered  matrix  S'  are  defined  as 


1  Ns 

Smj  =  Smj  -  VT 

5 1=1 


(7a) 


Similarly, 


(7b) 

5  1  =  1 

PCA  looks  for  a  matrix  P  that  decomposes  X  into  virtual 
sources,  S  =  PX.  It  also  holds  that  S'  =  PX ',  since  P  is  just  a 
rotation  matrix  which  does  not  change  the  center  of  the  data. 


cov(S)  =  S' S' T  =  (PX')(PX'f  =  PX'X'TPT  =  A,  (8) 


where  A  =  diagfAi,  A2, . . .}.  The  eigenvalues  Am  are  variances 
in  the  covariance  matrix.  Therefore,  X'X'TPT  =  PT A, 
where  PT  is  orthonormal.  PCA  is  realized  by  eigenvalue 
decomposition  (EVD)  of  the  covariance  matrix  of  X.  The 
eigenvectors  with  leading  eigenvalues  (largest  variances)  are 
selected  to  be  the  PCs  using  the  L-curve  [24]. 

Since  X  =  PTS  ~  AS ,  A  is  determined  as  a  matrix 
including  only  PCs.  S  is  calculated  as  S  «  ( AT A )  ATX. 
Rows  of  S  and  columns  of  A  represent  principal  component 
intensity  distributions  (PCIDs)  on  the  source  plane  and 
detector  plane,  respectively  and  are  proportional  to  the 
images  of  the  virtual  sources  projected  on  the  source  and 
detector  planes.  The  target  positions  are  determined  using 
(6a)  and  (6b). 


2.3.3.  NMF.  NMF  is  a  group  of  multivariate  analysis  algo¬ 
rithms  that  factorize  a  matrix  X  into  A,  and  S  :  X  =  AS, 
A  and  S  are  nonnegative  [6].  Unlike  ICA  and  PCA,  NMF 
does  not  imply  any  relationship  between  the  retrieved  com¬ 
ponents;  instead,  it  just  enforces  non-negativity  constraints 
on  A  and  S.  There  are  various  algorithms  developed  to  solve 
NMF,  such  as  the  multiplicative  update  method  [5]  and 
alternating  least  squares  method  [25,  26]. 

In  the  multiplicative  update  implementation  of  NMF,  A 
and  S  can  be  found  by  minimizing  the  square  of  Euclidean 
distance  ||X  -  AS||2  as  the  cost  function,  where  A  >  0  and 
S  >  0,  using  the  multiplicative  update  rule: 


Ah  < —  Aik 


C ASST)ik ’ 


(9a) 


Skj  ‘ —  Skj 


(ATX)kj 
C ATAS)kj ■ 


(9b) 


The  alternating  least  squares  implementation  of  NMF 
uses  alternate  least  squares  steps  to  estimate  A  (or  S), 
and  use  that  estimate  to  optimize  S  (or  A),  repeating  the 
alternative  steps  until  the  desired  optimization  is  obtained. 
Nonnegativity  is  ensured  by  setting  any  negative  element  of 
A  or  S  equal  to  0. 

An  NMF  toolbox  was  obtained  from  http://cogsys.imm 
.dtu.dk/toolbox/  to  perform  NMF  computation.  A  built-in 
command  nnmf  is  also  available  in  Matlab  (R201  la). 

NMF  algorithm  requires  that  the  non -negativity  assump¬ 
tion  must  hold  in  the  problem.  In  particular,  for  absorptive 
targets,  when  X  is  constructed  with  -A 0,  rm  should  be 
positive,  that  is,  the  targets  should  be  more  absorbing  than 
the  background.  If  the  targets  have  weaker  attenuation 
properties  than  the  background,  X  needs  to  be  constructed 
with  +A0  instead.  For  scattering  targets,  X  should  be  treated 
similarly  to  keep  its  elements  positive. 

When  NMF  is  applied  to  a  scattering  target,  only  the 
centrosymmetric  component  shows  up  properly,  since  the 
other  two  components  have  dumb-bell  shape  which  includes 
negative  values  [8].  So  without  any  prior  knowledge  or  some 
other  experimental  means  to  assess  if  the  target  is  absorptive 
or  scattering,  NMF  may  not  distinguish  between  the  two 
possibilities. 
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Figure  1:  Light  intensity  distribution  on  the  detector  plane  is 
recorded  when  a  point  source  scans  on  the  source  plane. 

The  decomposition  methods  can  be  applied  with  differ¬ 
ent  sample  geometries  such  as  slab  and  cylindrical  geome¬ 
tries,  and  different  measurement  domains  such  as  time- 
resolved  domain,  frequency  domain,  and  continuous  wave 
(CW).  In  this  paper,  Greens  functions  for  slab  geometry 
[23]  with  CW  measurement  were  used  for  simulation  and 
experiments. 

3.  Simulation 

The  sample  was  considered  to  be  a  40  mm  thick  uniform 
scattering  slab  with  lateral  dimension  of  80  mm  X  80  mm,  as 
shown  in  Figure  1.  Its  absorption  and  diffusion  coefficients 
were  taken  to  be  \ia  =  0.003  mm-1  and  D  =  1/3  mm  (trans¬ 
port  mean  free  path,  lt  =  1  mm),  respectively,  which  are 
similar  to  the  average  value  of  those  parameters  for  human 
breast  tissue.  An  absorptive  and  a  scattering  point  targets 
were  placed  at  (50,  60,  15)  mm  and  (30,  30,  25)  mm,  respec¬ 
tively.  The  index  of  refraction  n  of  the  medium  was  taken  to 
be  1.33.  The  speed  of  light  is  2.998  X  108  m/s  or  299.8  mm/ns 
in  vacuum,  and  225.4  mm/ns  in  the  medium.  The  absorption 
coefficient  of  the  absorptive  target  was  set  to  be  higher  than 
the  background  with  A[ia  =  0.001  mm-1,  while  the  diffusion 
coefficient  was  taken  to  be  the  same  as  that  of  background. 
The  diffusion  coefficient  of  the  scattering  target  was  set  to 
be  lower  than  the  background  (higher  scattering  coefficient) 
with  A D  =  -0.1mm  ( lt  =  0.7  mm),  while  the  absorption 
coefficient  was  taken  to  be  the  same  as  that  of  the  back¬ 
ground.  The  volumes  of  both  targets  are  set  to  be  8  mm3.  The 
optical  strengths  of  the  absorptive  and  scattering  targets  were 
AftacAV  =  1.803  mm3 /ns  and  ADc AV  =  -180.3  mm5/ns, 
respectively.  The  incident  CW  beam  step  scanned  the  sample 
at  21  X  21  grid  points  covering  an  80  X  80  mm2  area,  with 
a  step  size  of  4  mm.  Light  on  the  opposite  side  was  recorded 
at  41  X  41  grid  points  covering  the  same  area.  Multiplicative 
Gaussian  noise  of  5%  was  added  to  the  simulated  data.  The 
data  matrix  X  was  then  obtained  using  (2)  and  (3)  directly 
and  analyzed  using  the  three  different  algorithms. 

3.1.  ICA  Analysis.  One  independent  component  for  the 
absorptive  target  and  three  independent  components  for  the 
scattering  target  were  retrieved  by  ICA.  The  independent 
component  intensity  distributions  (ICIDs)  on  the  detector 


plane  are  shown  in  Figures  2(a),  2(c),  2(d),  and  2(e).  Similar 
ICIDs  were  obtained  on  the  source  plane.  Figure  2(g)  shows 
the  centrosymmetric  ICID  for  the  scattering  target,  and 
Figure  2(i)  shows  the  ICID  for  the  absorptive  target,  on  the 
source  plane. 

The  components  in  either  the  detector  plane  or  the 
source  plane  can,  in  principle,  be  used  to  extract  position 
and  optical  strength  of  the  target(s).  FFowever,  in  our 
experimental  arrangement  signal  is  collected  by  a  1024  X 
1024  pixels  CCD  camera,  while  the  source  plane  is  scanned 
in  an  x-y  array  of  points,  which  is  much  smaller  than  the 
number  of  pixels  in  the  CCD  camera.  Consequently,  the 
resolution  in  the  detector  plane  is  much  better,  and  the  data 
set  more  robust  than  the  source  side.  So,  we  used  the  images 
on  the  detector  plane  for  retrieving  target  information  using 
experimental  data.  While  it  would  not  matter  in  simulation, 
to  be  consistent  with  experimental  situations,  we  employed 
detector  plane  images  when  using  simulated  data  as  well  for 
all  three  algorithms.  Table  1  lists  the  locations  and  strengths 
of  the  absorptive  and  scattering  targets  retrieved  by  fitting  the 
spatial  intensity  profile  of  the  centrosymmetric  components 
on  the  detector  plane  to  Greens  functions  and  derivatives  of 
Greens  functions  using  (6a)  and  (6b),  respectively,  as  shown 
in  Figures  2(b)  and  2(f).  Figures  2(h)  and  2(j)  show  the 
corresponding  fits  to  the  profiles  on  the  source  plane. 

3.2.  PCA  Analysis.  Eigenvalue  equation  of  the  covariance 
matrix  of  X  was  solved.  The  eigenvalues  found  by  PCA  were 
sorted  in  a  descending  order.  Figure  3  shows  a  plot  of  leading 
20  eigenvalues  on  a  logarithmic  scale. 

First  four  leading  eigenvalues  were  selected  for  PCs. 
The  corresponding  PCIDs  were  calculated.  The  PCIDs 
on  the  detector  plane  are  shown  in  Figure  4.  Similar 
images  for  PCIDs  on  the  source  plane  were  obtained.  The 
scattering  target  has  one  centrosymmetric  (Figure  4(a)) 
component  and  two  dumb-bell  shaped  (Figures  4(c)  and 
4(d))  components,  while  the  absorptive  target  has  only  one 
component  (Figure  4(e)). 

Figures  4(b)  and  4(f)  show  fits  to  the  spatial  intensity 
profile  of  the  centrosymmetric  component  of  the  scattering 
target  and  that  of  the  absorptive  target,  respectively,  to 
retrieve  the  locations  of  the  two  targets.  The  locations  and 
optical  strengths  of  the  targets  retrieved  by  PCA  are  also 
shown  in  Table  1. 

3.3.  NMF  Analysis.  The  mixing  matrix  and  virtual  sources 
were  retrieved  from  the  data  matrix  X  using  NMF  as 
explained  in  Section  2.3.3.  As  in  the  other  two  approaches, 
only  one  component  is  extracted  for  the  absorptive  tar¬ 
get.  Since  NMF  has  a  non -negativity  constraint,  only  the 
centrosymmetric  component  for  the  scattering  target  is 
obtained.  Nonnegative  component  intensity  distributions 
(NCIDs)  on  detector  planes  are  shown  in  Figure  5.  Similar 
images  for  NCIDs  on  source  plane  were  also  obtained  using 
the  virtual  sources  in  S.  The  results  are  also  shown  in  Table  1. 

3.4.  Results  and  Discussion.  The  positions  and  optical 
strengths  of  the  targets  retrieved  by  ICA,  PCA,  and  NMF 
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Figure  2:  ICA-extracted  two-dimensional  intensity  distribution  on  the  detector  plane  of:  (a)  the  centrosymmetric  component;  (c)  and  (d) 
dumb-bell  shaped  components  of  the  scattering  target;  (e)  the  absorptive  target.  Similar  intensity  distribution  on  the  source  plane  of:  (g) 
the  centrosymmetric  component  of  the  scattering  target  and  (i)  the  absorptive  target  for  comparison.  Fits  to  the  spatial  intensity  profile  on 
the  detector  plane  along  the  white  dashed  line  (shown  in  figures)  of  the  centrosymmetric  component  of  the  scattering  target  is  shown  in 
(b),  and  that  of  the  absorptive  target  is  shown  in  (f).  Corresponding  fits  to  spatial  profiles  on  the  source  plane  are  displayed  in  (h)  and  (j), 
respectively. 


Figure  3:  A  logarithmic  plot  of  the  first  20  PCA  eigenvalues. 


algorithms  are  shown  in  Table  1,  and  compared  to  the  known 
values.  The  retrieved  results  using  all  three  algorithms  from 
this  simulated  data  are  in  excellent  agreement  with  the 
known  values. 


4.  Experiments 

4.1.  Experimental  Materials  and  Methods.  In  this  Section, 
the  algorithms  are  evaluated  using  experimental  data  for 
absorptive  and  scattering  targets  embedded  in  model  scat¬ 
tering  media  whose  absorption  and  scattering  properties  are 
adjusted  to  mimic  the  average  values  of  those  parameters 
for  human  breast  tissues.  Two  different  experiments  were 
carried  out  with  two  different  samples.  The  first  sample 
used  a  250  mm  X  250  mm  X  50  mm  transparent  plastic 
container  filled  with  Intralipid  10%  suspension  in  water  as 
the  background  medium.  The  concentration  of  Intralipid 
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Table  1:  Positions  and  optical  strengths  retrieved  using  simulated  data  and  ICA,  PCA,  and  NMF  algorithms. 


Target 

Known  position  (mm) 

Algorithm 

Fitted  position  (mm) 

Error  (mm) 

Known  strength* 

Fitted  strength* 

Error  (%) 

ICA 

(29.9,  30.0,25.1) 

(0.1,0,  0.1) 

-180.3 

-179.9 

0.22 

Sea. 

(30,  30,  25) 

PCA 

(30.0,  30.0,  25.0) 

(0,  0,  0) 

-180.3 

-180.1 

0.11 

NMF 

(30.0,  30.0,  25.0) 

(0,  0,  0) 

-180.3 

-178.5 

1 

ICA 

(50.1,60.2,  15.0) 

(0.1,  0.2,  0) 

1.803 

1.826 

1.28 

Abs. 

(50,  60,  15) 

PCA 

(50.1,60.1,  14.9) 

(0.1,  0.1,  0.1) 

1.803 

1.812 

0.5 

NMF 

(50.1,60.1,  15.0) 

(0.1,  0.1,0) 

1.803 

1.803 

0 

*The  unit  for  absorption  strength  of  the  target  is  mm3 /ns  and  for  scattering  strength  is  mm5 /ns. 


10%  was  adjusted  to  provide  [27,  28]  an  absorption  coef¬ 
ficient  of  [ia  ~  0.003  mm-1,  and  a  transport  mean-free 
path  lt  ~  1.43  mm  at  785  nm.  The  second  sample  used  a 
similar  container  with  dimension  of  250  mm  X  250  mm  X 
60  mm  filled  with  Intralipid  20%  suspension  in  water.  The 
concentration  of  Intralipid  20%  was  adjusted  to  provide 
[27,  28]  [ ia  ~  0.003  mm-1,  and  lt  ~  1  mm  at  785  nm.  These 
optical  parameters  of  the  medium  were  selected  to  be  similar 
to  the  average  values  of  those  parameters  for  human  breast 
tissue.  The  thickness  of  the  samples  was  also  comparable  to 
that  of  a  typical  compressed  female  human  breast. 

In  the  first  experiment,  two  absorptive  targets  were 
embedded  in  the  medium.  The  targets  were  ~10-mm 
diameter  glass  spheres  filled  Indocyanine  green  (ICG)  dye 
dissolved  in  Intralipid-20%  suspension  in  water  to  obtain  an 
absorption  coefficient  [ ia  =  1.15  mm-1  at  785  nm,  and  to 
match  the  background  scattering  coefficient  of  1.97  mm-1. 
The  targets  were  placed  at  (57.2,  18.1,  20.0)  mm  and  (19.9, 
48.1,  25.0)  mm,  respectively. 

In  the  second  experiment,  two  scattering  targets  were 
embedded,  which  were  also  ~10mm  diameter  glass  spheres, 
filled  with  Intralipid-20%  suspension  in  water.  The  transport 
mean  free  path,  lt  was  adjusted  to  be  0.25  mm,  with  scattering 
coefficient  [is  «  1 1  mm-1,  and  absorption  coefficient  [ia  same 
as  that  of  the  background  medium.  The  targets  were  placed 
in  the  middle  plane  (z  =  30  mm)  in  the  container  with  a 
lateral  distance  of  40  mm  from  each  other  (center  to  center). 

The  experimental  setup  is  schematically  shown  in 
Figure  6.  A  10  mW  785  nm  diode  laser  beam  was  used 
to  illuminate  the  first  sample,  while  a  100  mW  785  nm 
diode  laser  beam  was  used  for  the  second  sample.  The 
input  surface  (source  plane)  of  the  samples  was  scanned 
across  the  laser  beam  in  an  x-y  array  of  grid  points  to 
realize  the  multi- source  interrogation  of  the  samples.  The 
transmitted  light  from  the  exit  surface  (detector  plane) 
was  recorded  by  a  1024  pixel  x  1024  pixel  (pixel  size  = 
24  fim)  CCD  camera  (Photometries  CFF350)  equipped  with 
a  60  mm  focal-length  camera  lens.  Each  pixel  of  the  CCD 
camera  can  be  considered  to  be  a  detector  implementing 
the  multidetector  signal  acquisition  arrangement.  A  set  of 
16  bit  1024  pixel  X  1024  pixel  images  were  acquired.  The  two 
samples  were  scanned  in  an  array  of  11  X  12  and  11  X  15 
grid  points,  respectively,  with  a  step  size  of  5  mm  in  both 
cases.  The  processes  of  scanning  and  data  acquisition  were 
controlled  by  a  personal  computer.  At  all  scan  positions,  raw 
transillumination  images  of  the  samples  were  recorded  by  the 
computer  for  further  analysis. 


4.2.  Analysis  and  Results.  A  region  of  interest  (ROI)  was 
cropped  out  from  each  image.  Then,  every  5x5  pixels  in 
each  cropped  image  were  binned  to  one  pixel  to  enhance 
signal- to -noise  ratio.  A  background  image  was  generated 
by  calculating  an  average  image  for  all  scan  positions  to 
approximate  the  transillumination  image  without  target(s) 
embedded. 

This  averaging  method  for  generating  background  image 
is  suitable  for  small  targets  used  in  our  experiments,  as 
the  ratio  of  the  volume  of  the  sample  to  that  of  the 
target  was  quite  high  (~500:1).  For  in  vivo  imaging  of 
tumors  in  early  stages  of  growth,  the  breast- to -tumor  volume 
ratio  will  be  similarly  high,  and  the  averaging  method 
will  be  applicable.  Alternative  approaches  for  generating  a 
background  image  include  using  image  of  (a)  a  phantom 
that  has  the  same  average  optical  properties  as  the  sample 

[29] ;  (b)  the  healthy  contralateral  breast  for  breast  imaging 

[30] ;  (c)  the  sample  obtained  using  light  of  wavelength  for 
which  the  target (s)  and  the  background  have  identical  optical 
properties  [31].  Still  another  approach  is  to  compute  the 
background  using  an  appropriate  forward  model  [32].  A 
more  detailed  discussion  of  this  important  issue  appears  in 
one  of  our  earlier  publications  [33]. 

The  background  image  was  also  cropped  and  binned  cor¬ 
responding  to  the  ROI  for  each  scan  position.  Perturbation 
in  the  light  intensity  distribution,  A 0  due  to  targets  in  each 
image  was  found  by  subtracting  the  background  image  from 
the  image.  The  data  matrix  X  was  then  constructed  using  the 
light  intensity  perturbations  at  all  scan  positions.  ICA,  PCA, 
and  NMF  decomposition  algorithms  were  performed  on 
the  data  matrix  separately.  Results  are  shown  and  discussed 
below. 

4.2.1.  Absorptive  Targets.  The  images  on  the  detector  plane 
obtained  using  the  ICA,  PCA,  and  NMF  algorithms  are 
shown  in  Figures  7,  8,  and  9,  respectively.  Similar  images 
on  the  source  plane  were  also  obtained  using  all  three  algo¬ 
rithms.  The  right  side  of  each  figure  shows  the  corresponding 
spatial  intensity  profile.  Focations  of  the  targets  are  extracted 
from  fits  to  these  spatial  intensity  profiles,  as  described  in 
Section  2.3  using  (6a)  and  (6b).  The  results  are  presented  in 
Table  2.  In  Figure  7,  images  on  the  source  plane  are  shown  in 
(e)  and  (g),  and  Greens  function  fits  to  their  spatial  profiles 
are  shown  in  (f)  and  (h)  for  comparison. 

It  follows  from  the  comparison  of  the  results  in  Table  2 
that  the  positions  retrieved  by  all  three  algorithms  are 
in  good  agreement  with  the  known  positions.  The  errors 
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Figure  4:  PCA-extracted  two-dimensional  intensity  distribution 
on  the  detector  plane  of:  (a)  the  centrosymmetric  component; 
and  (c)  and  (d)  dumb-bell  shaped  components  of  the  scattering 
target;  (e)  the  absorptive  target.  Green  s  function  fits  to  the  spatial 
intensity  profiles  along  the  dashed  line  (shown  in  figures)  of  the 
(b)  centrosymmetric  component  of  the  scattering  target  and  (f) 
absorptive  target,  respectively,  to  retrieve  the  locations  of  the  two 
targets. 


in  the  retrieved  locations  (x,y,z)  of  the  two  targets  were 
within  1.7  mm.  The  PCIDs  were  not  totally  separated.  Some 
“residue”  was  observed  in  one  PCID  from  the  other.  ICA 
and  NMF  separated  two  components  from  this  dataset  more 
clearly. 


4.2.2.  Scattering  Targets.  The  “images”  corresponding  to 
the  centrosymmetric  components  of  the  virtual  sources 
(targets)  on  the  detector  plane  obtained  using  the  ICA, 
PCA,  and  NMF  algorithms  are  shown  in  Figures  10,  11, 
and  12,  respectively.  Similar  images  on  the  source  plane 
were  also  obtained.  The  right  side  of  each  figure  shows 
the  corresponding  spatial  intensity  profile.  Locations  of  the 
targets  are  extracted  from  fits  to  these  spatial  intensity 


Figure  5:  NMF-extracted  two-dimensional  intensity  distribution 
on  the  detector  plane  of:  (a)  the  centrosymmetric  component  of  the 
scattering  target;  (c)  the  absorptive  target.  Fits  to  the  corresponding 
spatial  intensity  profiles  along  the  dashed  line  (shown  in  figures)  are 
shown  in  (b)  and  (d),  respectively. 


Figure  6:  A  schematic  diagram  of  the  experimental  arrangement 
used  for  imaging  objects  embedded  in  a  turbid  medium.  The  inset 
at  the  bottom  shows  the  2D  array  in  the  input  plane  that  was 
scanned  across  the  incident  laser  beam;  the  inset  to  the  right  shows 
a  typical  raw  image  recorded  by  the  CCD.  (CCD:  charge  coupled 
device,  and  PC:  personal  computer). 


profiles,  as  described  in  Section  2.3  using  (6a)  and  (6b).  The 
results  are  presented  in  Table  3. 

Both  targets  were  detected  by  all  three  algorithms.  The 
target  locations  retrieved  by  three  algorithms  are  shown 
in  Table  3  and  compared  with  known  locations.  Overall, 
all  three  algorithms  detect  and  locate  the  scattering  and 
the  absorptive  targets  with  good  accuracy,  the  maximum 
deviation  of  any  one  coordinate  from  the  known  value 
being  ~3mm.  Since  the  maximum  difference  between  the 
known  and  retrieved  position  coordinates  was  larger  for 
the  scattering  targets,  we  calculated  the  squared  correlation 
coefficient  y  to  assess  the  fitting  quality.  NMF  retrieves 
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Detector  plane  images  and  profiles 
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Figure  7:  ICA-generated  ICIDs  on  the  detector  plane  are  shown  in 
(a)  and  (c);  corresponding  Greens  function  fits  to  the  horizontal 
spatial  profiles  through  the  dashed  lines  are  shown  in  (b)  and  (d). 
ICIDs  on  the  source  plane  are  shown  in  (e)  and  (g);  corresponding 
Green’s  function  fits  to  the  horizontal  spatial  profiles  through  the 
dashed  line  are  shown  in  (f)  and  (h). 


the  position  coordinates  better  (within  0.5  mm)  for  the  scat¬ 
tering  Target  2  than  done  by  ICA  and  PCA  (deviation  from 
known  values  being  between  2-3  mm).  NMF  retrieved  the 
position  coordinates  for  Target  1  with  3.0  mm  error  in  z 
direction,  which  is  not  as  good  as  that  done  by  ICA  and 
PCA.  But  y  is  0.783  and  0.778  in  the  fittings  for  ICA 
and  PCA,  respectively,  as  compared  to  0.993  for  NMF, 
indicating  that  the  quality  of  the  fitting  is  better  for  NMF. 
The  quality  of  fitting  is  presumably  affected  by  the  efficacy  of 
decomposition.  The  decomposed  NCIDs  by  NMF  were  more 
“clean”  than  those  decomposed  by  ICA  and  PCA.  We  ascribe 
the  observed  higher  errors  in  ICA  and  PCA  estimates  of  the 


xlO4 
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x  (mm) 


(c) 


(d) 


Figure  8:  PCIDs  on  the  detector  plane  are  shown  in  (a)  and  (c);  and 
corresponding  Green’s  function  fits  to  the  horizontal  spatial  profiles 
through  the  dashed  line  are  shown  in  (b)  and  (d). 


Figure  9:  NCIDs  on  the  detector  plane  are  shown  in  (a)  and  (c); 
corresponding  Green’s  function  fits  to  the  horizontal  spatial  profiles 
through  the  dashed  line  are  shown  in  (b)  and  (d). 


position  coordinates  of  the  scattering  Target  2  than  the  NMF 
estimates  to  the  interference  from  the  other  virtual  source 
(corresponding  to  Target  1)  in  ICA  (Figure  10(c))  and  PCA 
(Figure  11(c))  images.  It  is  commonly  believed  that  errors  in 
locating  a  scattering  target  are  higher  than  that  for  locating 
an  absorptive  target,  and  the  results  of  this  study  conform  to 
that  notion. 
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Table  2:  Known  positions  versus  retrieved  positions  of  the  absorptive  targets  using  ICA,  PCA,  and  NMF  algorithms. 


Target 

Known  position  (mm) 

Algorithm 

Fitted  position  (mm) 

Error  (mm) 

ICA 

(57.4,  18.2,21.5) 

(0.2,  0.1,  1.5) 

1 

(57.2,  18.1,20) 

PCA 

(57.4,  18.2,  20.6) 

(0.2,  0.1,  0.6) 

NMF 

(57.4,  18.2,  19.5) 

(0.2,  0.1,  0.5) 

ICA 

(18.2,  46.7,  24.7) 

(1.7,  1.4,  0.3) 

2 

(19.9,  48.1,25) 

PCA 

(18.2,  47.6,25.9) 

(1.7,  0.5,  0.9) 

NMF 

(18.2,  47.6,23.3) 

(1.7,  0.5,  1.7) 

Table  3:  Known  positions  versus  retrieved  positions  of  the  scattering  targets  using  ICA,  PCA,  and  NMF  algorithms. 


Target 

Known  position  (mm) 

Algorithm 

Fitted  position  (mm) 

Error  (mm) 

ICA 

(12.6,28.7,  29.1) 

(0.4,  0.7,  0.9) 

1 

(13.0,  28.0,  30.0) 

PCA 

(12.6,28.7,28.6) 

(0.4,07,  1.4) 

NMF 

(12.0,28.5,  33.0) 

(1.0,  0.5,  3.0) 

ICA 

(51.0,31.8,  26.8) 

(2.3,  3.3,  3.2) 

2 

(53.3,28.5,30.0) 

PCA 

(50.9,31.8,  26.7) 

(2.4,  3.3,  3.3) 

NMF 

(53.3,  28.0,  30.3) 

(0.0,  0.5,  0.3) 
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Figure  10:  ICA-generated  ICIDs  on  the  detector  plane  are  shown 
in  (a)  and  (c);  corresponding  Green  s  function  fits  to  the  horizontal 
spatial  profiles  through  the  dashed  line  are  shown  in  (b)  and  (d). 


Figure  11:  PCIDs  on  the  detector  plane  are  shown  in  (a)  and  (c); 
corresponding  Greens  function  fits  to  the  horizontal  spatial  profiles 
through  the  dashed  line  are  shown  in  (b)  and  (d). 


5.  Summary  and  Discussion 

Diffusive  optical  imaging  was  modeled  as  a  BSS  problem. 
ICA,  PCA,  and  NMF  were  used  to  decompose  the  data  matrix 
and  locate  the  targets  embedded  in  a  highly  scattering  turbid 
medium.  Only  the  components  corresponding  to  the  targets 
were  extracted  from  a  large  dataset  for  target  detection  and 
localization. 

It  may  be  instructive  to  compare  the  objectives, 
scope,  and  computational  complexity  of  these  decompo¬ 
sition  methods  with  model-based  reconstruction  methods. 
Decomposition  methods  obtain  the  3D  locations  of  tar¬ 
gets  (the  number  of  targets  is  generally  small).  Based  on 
the  retrieved  locations,  the  methods  may  then  be  further 


extended  to  retrieve  size  and  optical  property  information 
of  the  targets  [9].  The  common  practice  of  model-based 
inverse  reconstruction  methods  is  to  discretize  the  sample 
volume  into  NxNxN  voxels  and  estimate  absorption  and/or 
scattering  coefficient  in  each  voxel  iteratively.  Voxels  with 
significantly  different  optical  properties  than  the  surround¬ 
ing  are  regions  of  interest  and  may  be  identified  as  targets. 
While  estimating  the  optical  properties,  the  forward  model 
is  solved  repeatedly  to  calculate  the  intensity  of  the  multiply 
scattered  light  on  the  sample  boundary.  The  difference 
between  the  intensity  of  the  multiply  scattered  light  predicted 
by  the  forward  model  and  the  experimental  measurements  is 
minimized  by  seeking  an  optimal  set  of  the  optical  properties 
of  every  voxel  in  the  sample  volume.  The  number  of  variables 
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Figure  12:  NCIDs  on  the  detector  plane  are  shown  in  (a)  and  (c); 
corresponding  Green  s  function  fits  to  the  horizontal  spatial  profiles 
through  the  dashed  line  are  shown  in  (b)  and  (d). 


is  thus,  on  the  order  of  N 3.  To  determine  location(s)  of 
target(s)  in  three  dimensions,  the  decomposition  methods 
process  the  data  matrix  to  retrieve  the  main  components  (A 
and  S).  Here,  A  and  S  are  two-dimensional  matrices  with  the 
number  of  unknowns  on  the  order  of  N2.  The  number  of 
unknowns  is,  hence,  reduced  N  times  in  the  decomposition 
methods  compared  to  the  model-based  approaches,  which 
leads  to  a  substantial  saving  in  the  computational  time 
when  N  is  large.  No  repeated  solution  of  the  forward 
model  is  involved  in  decomposition  methods.  Consequently, 
decomposition  methods  are  considerably  faster. 

A  comparison  of  the  computational  complexity  of  these 
two  types  of  approaches  may  shed  further  light  on  their 
relative  computation  economy.  For  a  model-based  iterative 
reconstruction  method,  an  equation  of  the  form  b  =  Wx 
is  solved  to  find  the  targets,  where  W  is  a  weight  matrix 
of  size  NdNs  x  Nv,  Nd,  Ns,  and  Nv  are  the  numbers  of 
detectors,  sources  and  voxels,  respectively,  b  is  an  NdNs  x 
1  vector  describing  the  perturbation  in  the  detected  light 
intensity  due  to  the  presence  of  targets,  and  x  is  the 
perturbation  in  the  optical  properties  from  the  background 
values  with  dimension  of  Nv  X  1.  The  computational 
complexity  is  typically  0(NdNsN2)  for  a  single  iteration.  For 
the  decomposition  approach,  b  is  written  as  a  2D  matrix 
X  with  dimension  Nd  X  Ns.  To  decompose  matrix  X,  the 
computational  complexity  per  iteration  is  typically  of  order 
0(NdNk )  for  ICA  [34],  and  0(NdNsNk )  for  NMF  [16],  where 
Nk  is  the  number  of  components  that  relates  to  the  number 
of  targets  and  is  usually  a  small  number.  For  PCA  using 
SVD,  the  complexity  is  0(N2Nk )  [34].  The  computational 
complexity  of  the  intrinsic  iterative  process  involved  in  the 
matrix  decomposition  algorithms  is  much  lower  than  that  in 
the  model-based  inverse  reconstruction  methods. 

All  three  matrix  decomposition  methods  presented  in 
this  manuscript  can  potentially  be  used  in  in-vivo  real¬ 
time  breast  cancer  imaging.  The  three  algorithms  have 


different  assumptions,  which  may  lead  to  different  favored 
conditions.  In  this  study,  the  algorithms  were  evaluated  using 
simulative  and  experimental  data  using  model  scattering 
media  and  absorptive  and  scattering  targets.  The  (x,y,z) 
positions  of  the  targets  were  retrieved  with  good  accuracy. 
The  decomposition  provided  by  ICA  is  “cleaner”  than  that 
of  the  PCA.  PCA  did  not  clearly  separate  the  two  absorptive 
targets  used  in  the  first  experiment.  NMF  decomposition 
seems  to  provide  residue-free  “cleaner”  images  than  the  other 
two  methods  in  this  study.  However,  since  NMF  is  based 
on  nonnegativity  assumption,  the  results  might  deteriorate 
when  such  a  non-negativity  assumption  does  not  hold  well. 
While  continuous  wave  measurements  were  used  in  the  work 
presented  in  this  paper,  the  approaches  could  be  used  with 
frequency  domain  and  time  domain  measurements  as  well. 

The  work  presented  here  focuses  on  detecting  and 
locating  small  targets,  which  derive  impetus  from  the  need  to 
detect  tumors  in  early  stages  of  growth  when  those  are  more 
amenable  to  treatment.  All  three  methods  are  applicable  for 
extended  targets  as  well  and  are  expected  to  provide  the 
“center  of  optical  strength”  as  the  location  of  the  target. 

All  three  approaches  are  applicable  for  both  scattering 
and  absorbing  targets  and  may  be  used  in  clinical  setting. 
The  contrast  between  a  tumor  and  surrounding  normal 
tissue  can  be  due  to  differences  in  absorption,  scattering,  or 
both  absorption  and  scattering  properties  and  may  depend 
significantly  on  the  wavelength  of  light  used.  However,  a 
priori  knowledge  of  the  optical  characteristics  (absorptive 
or  scattering)  is  not  crucial.  As  has  been  shown  in  (2) 
and  (3),  the  expression  for  elements  of  the  data  matrix  for 
absorptive  targets  involves  Greens  Functions  G,  while  that 
for  scattering  targets  involves  dG/dz  «  - kG ,  where  k  = 
^pJD  in  CW  [9].  This  relationship  with  G  provides  basis 
for  detection  and  localization  of  target(s),  whether  contrast 
is  due  to  absorption,  scattering,  or  both.  We  are  using 
transillumination  geometry,  which  is  one  of  the  approaches 
used  by  other  researchers,  and  adequate  signal  for  in  vivo 
breast  imaging  is  obtained  [29,  35-38]. 

In  this  paper,  we  presented  results  when  the  approaches 
were  used  to  detect  and  obtain  three-dimensional  location 
information  of  the  targets.  We  have  demonstrated,  while 
developing  OPTICA  that  a  backprojection  formalism  can 
be  further  implemented  to  get  a  cross-section  image  of  the 
target  [11],  or  the  retrieved  target  locations  can  be  fed  into 
other  DOI  methods  as  a  priori  information  to  get  three- 
dimensional  tomographic  images.  Since  the  approaches  are 
suited  for  small  targets,  this  hold  promise  for  detecting  and 
locating  breast  tumors  in  early  stages  of  growth,  which  is 
crucially  important  for  effective  treatment.  Further  work 
involving  ex  vivo  (model)  and  in  vivo  imaging  of  cancerous 
breast  will  be  needed  to  establish  the  full  potential  of  these 
approaches. 


Acknowledgment 

The  paper  is  supported  in  part  by  US  Army  Medical  Research 
and  Materiel  Command  under  Contract  no.  W81XWH-07- 
1-0454. 


International  Journal  of  Optics 


Appendix  1 


11 


References 

[  1  ]  X.  R.  Cao  and  R.  W.  Liu,  “General  approach  to  blind  source 
separation,”  IEEE  Transactions  on  Signal  Processing ,  vol.  44,  no. 
3,  pp.  562-571,  1996. 

[2]  J.  R  Cardoso,  “Blind  signal  separation:  statistical  principles,” 
Proceedings  of  the  IEEE,  vol.  86,  no.  10,  pp.  2009-2025,  1998. 

[3]  A.  Hyvarinen,  J.  Karhunen,  and  E.  Oja,  Independent  Compo¬ 
nent  Analysis,  Wiley,  New  York,  NY,  USA,  2001. 

[4]  I.  T.  Jolliffe,  Principal  Component  Analysis,  Springer,  New  York, 
NY,  USA,  1986. 

[5]  D.  D.  Lee  and  H.  S.  Seung,  “Learning  the  parts  of  objects  by 
non-negative  matrix  factorization,”  Nature,  vol.  401,  no.  6755, 
pp.  788-791,  1999. 

[6]  M.  W.  Berry,  M.  Browne,  A.  N.  Langville,  V.  R  Pauca,  and  R. 
J.  Plemmons,  “Algorithms  and  applications  for  approximate 
nonnegative  matrix  factorization,”  Computational  Statistics 
and  Data  Analysis,  vol.  52,  no.  1,  pp.  155-173,  2007. 

[7]  M.  Xu,  M.  Alrubaiee,  S.  K.  Gayen,  and  R.  R.  Alfano,  “Optical 
imaging  of  turbid  media  using  independent  component 
analysis:  theory  and  simulation,”  Journal  of  Biomedical  Optics, 
vol.  10,  no.  5,  Article  ID  051705,  2005. 

[8]  M.  Xu,  M.  Alrubaiee,  S.  K.  Gayen,  and  R.  R.  Alfano,  “Three- 
dimensional  localization  and  optical  imaging  of  objects  in 
turbid  media  with  independent  component  analysis,”  Applied 
Optics,  vol.  44,  no.  10,  pp.  1889-1897,  2005. 

[9]  M.  Xu,  M.  Alrubaiee,  S.  K.  Gayen,  and  R.  R.  Alfano,  “Optical 
diffuse  imaging  of  an  ex  vivo  model  cancerous  human  breast 
using  independent  component  analysis,”  IEEE  Journal  on 
Selected  Topics  in  Quantum  Electronics,  vol.  14,  no.  1,  pp.  43- 
49,  2008. 

[10]  M.  Alrubaiee,  M.  Xu,  S.  K.  Gayen,  M.  Brito,  and  R.  R. 
Alfano,  “Three-dimensional  optical  tomographic  imaging  of 
scattering  objects  in  tissue- simulating  turbid  media  using 
independent  component  analysis,”  Applied  Physics  Letters,  vol. 
87,  no.  19,  Article  ID  191112,  pp.  1-3,  2005. 

[11]  M.  Alrubaiee,  M.  Xu,  S.  K.  Gayen,  and  R.  R.  Alfano, 
“Localization  and  cross  section  reconstruction  of  fluorescent 
targets  in  ex  vivo  breast  tissue  using  independent  component 
analysis,”  Applied  Physics  Letters,  vol.  89,  no.  13,  Article  ID 
133902,  2006. 

[12]  D.  L.  Massart,  B.  G.  M.  Vandeginste,  S.  M.  Deming,  Y. 
Michotte,  and  L.  Kaufman,  Chemometrics:  A  Textbook,  Else¬ 
vier,  Amsterdam,  The  Netherlands,  1988. 

[13]  M.  Turk  and  A.  Pentland,  “Eigenfaces  for  recognition,”  Journal 
of  Cognitive  Neuroscience,  vol.  3,  no.  1,  pp.  71-86,  1991. 

[14]  L.  K.  Hansen,  J.  Larsen,  F.  A.  Nielsen  et  al.,  “Generalizable 
patterns  in  neuroimaging:  how  many  principal  components?” 
Neuroimage,  vol.  9,  no.  5,  pp.  534-544,  1999. 

[15]  J.  P.  Brunet,  P.  Tamayo,  T.  R.  Golub,  and  J.  P.  Mesirov, 
“Metagenes  and  molecular  pattern  discovery  using  matrix 
factorization,”  Proceedings  of  the  National  Academy  of  Sciences 
of  the  United  States  of  America,  vol.  101,  no.  12,  pp.  4164-4169, 
2004. 

[16]  V.  P.  Pauca,  J.  Piper,  and  R.  J.  Plemmons,  “Nonnegative  matrix 
factorization  for  spectral  data  analysis,”  Linear  Algebra  and  Its 
Applications,  vol.  416,  no.  1,  pp.  29-47,  2006. 

[17]  A.  Ishimaru,  “Diffusion  of  a  pulse  in  densely  distributed 
scatterers,”  Journal  of  the  Optical  Society  of  America,  vol.  68, 
no.  8,  pp.  1045-1050,  1978. 

[18]  K.  Furutsu,  “Diffusion  equation  derived  from  the  space- time 
transport  equation,”  Journal  of  the  Optical  Society  of  America, 
vol.  70,  no.  4,  pp.  360-366,  1980. 

[19]  M.  S.  Patterson,  B.  Chance,  and  B.  C.  Wilson,  “Time 
resolved  reflectance  and  transmittance  for  the  non-invasive 


measurement  of  tissue  optical  properties,”  Applied  Optics,  vol. 
28,  no.  12,  pp.  2331-2336,  1989. 

[20]  S.  Chandrasekhar,  Radiative  Transfer,  Clarendon  Press, 
Oxford,  UK,  1950. 

[21]  A.  Ishimaru,  Wave  Propagation  and  Scattering  in  Random 
Media,  vol.  1  of  Single  Scattering  and  Transport  Theory, 
Academic,  New  York,  NY,  USA,  1978. 

[22]  S.  R.  Arridge  and  J.  C.  Schotland,  “Optical  tomography: 
forward  and  inverse  problems,”  Inverse  Problems,  vol.  25, 
Article  ID  123010,  2009. 

[23]  S.  R.  Arridge,  “Photon-measurement  density  functions.  Part 
I:  analytical  forms,”  Applied  Optics,  vol.  34,  no.  31,  pp.  7395- 
7409,  1995. 

[24]  P.  C.  Hansen,  “Analysis  of  discrete  ill-posed  problems  by 
means  of  the  L-curve,”  SIAM  Review,  vol.  34,  pp.  561-580, 
1992. 

[25]  P.  Paatero  and  U.  Tapper,  “Positive  matrix  factorization:  a 
non-negative  factor  model  with  optimal  utilization  of  error 
estimates  of  data  values,”  Environmetrics,  vol.  5,  no.  2,  pp.  111- 
126,  1994. 

[26]  P.  Paatero,  “The  multilinear  engine:  a  table-driven,  least 
squares  program  for  solving  multilinear  problems,  including 
the  n-way  parallel  factor  analysis  model,”  Journal  of  Computa¬ 
tional  and  Graphical  Statistics,  vol.  8,  no.  4,  pp.  854-888,  1999. 

[27]  H.  J.  van  Staveren,  C.  J.  M.  Moes,  J.  van  Marie,  S.  A.  Prahl,  and 
M.  J.  C.  van  Gemert,  “Light  scattering  in  Intralipid- 10%  in  the 
wavelength  range  of  400-1100  nm,”  Applied  Optics,  vol.  3 1 ,  pp. 
4507-4514,  1991. 

[28]  C.  Bordier,  C.  Andraud,  E.  Charron,  J.  Lafait,  M.  Anastasi- 
adou,  and  A.  De  Martino,  “Illustration  of  a  bimodal  system  in 
Intralipid-20%  by  polarized  light  scattering:  experiments  and 
modeling,”  Applied  Physics  A,  vol.  94,  no.  2,  pp.  347-355, 2009. 

[29]  T.  Nielsen,  B.  Brendel,  R.  Ziegler  et  al.,  “Linear  image 
reconstruction  for  a  diffuse  optical  mammography  system  in 
a  noncompressed  geometry  using  scattering  fluid,”  Applied 
Optics,  vol.  48,  no.  10,  pp.  D1-D13,  2009. 

[30]  Y.  Ardeshirpour,  M.  Huang,  and  Q.  Zhu,  “Effect  of  the  chest 
wall  on  breast  lesion  reconstruction,”  Journal  of  Biomedical 
Optics,  vol.  14,  no.  4,  p.  044005,  2009. 

[31]  B.  W.  Pogue,  M.  S.  Patterson,  H.  Jiang,  and  K.  D.  Paulsen, 
“Initial  assessment  of  a  simple  system  for  frequency  domain 
diffuse  optical  tomography,”  Physics  in  Medicine  and  Biology, 
vol.  40,  no.  10,  pp.  1709-1729,  1995. 

[32]  A.  Poellinger,  J.  C.  Martin,  S.  L.  Ponder  et  al.,  “Near- infrared 
laser  computed  tomography  of  the  breast.  First  clinical 
experience,”  Academic  Radiology,  vol.  15,  no.  12,  pp.  1545- 
1553,2008. 

[33]  B.  Wu,  W.  Cai,  M.  Alrubaiee,  M.  Xu,  and  S.  K.  Gayen, 
“Time  reversal  optical  tomography:  locating  targets  in  a  highly 
scattering  turbid  medium,”  Optics  Express,  vol.  19,  no.  22,  pp. 
21956-21976,  2011. 

[34]  T.  Ristaniemi  and  J.  Joutsensalo,  “Advanced  ICA-based 
receivers  for  block  fading  DS-CDMA  channels,”  Signal  Process¬ 
ing,  vol.  82,  no.  3,  pp.  417-431,  2002. 

[35]  M.  Cutler,  “Transillumination  as  an  aid  in  the  diagnosis  of 
breast  lesions,”  Surgery,  Gynecology  and  Obstetrics,  vol.  48,  pp. 
721-729,  1929. 

[36]  R.  Choe,  A.  Corlu,  K.  Lee  et  al.,  “Diffuse  optical  tomography  of 
breast  cancer  during  neoadjuvant  chemotherapy:  a  case  study 
with  comparison  to  MRI ,”  Medical  Physics,  vol.  32,  no.  4,  pp. 
1128-1139,  2005. 

[37]  A.  Pifferi,  P.  Taroni,  A.  Torricelli,  F.  Messina,  R.  Cubeddu, 
and  G.  Danesini,  “Four- wavelength  time-resolved  optical 


12 


Appendix  1 


International  Journal  of  Optics 


mammography  in  the  680-980-nm  range,”  Optics  Letters ,  vol. 
28,  no.  13,  pp.  1138-1140,  2003. 

[38]  B.  J.  Tromberg,  B.  W.  Pogue,  K.  D.  Paulsen,  A.  G.  Yodh,  D.  A. 
Boas,  and  A.  E.  Cerussi,  “Assessing  the  future  of  diffuse  optical 
imaging  technologies  for  breast  cancer  management,”  Medical 
Physics ,  vol.  35,  no.  6,  pp.  2443-2451,  2008. 


Appendix  2 


Time-reversal  optical  tomography:  detecting  and  locating  extended 

targets  in  a  turbid  medium 

Binlin  Wu*ab,  W.  Caia,  M.  Xuc,  and  S.  K.  Gayenab 

aPhysics  Department,  The  City  College  of  the  City  University  of  New  York,  160  Convent  Ave,  New 

York,  NY  10031; 

bThe  Graduate  Center  of  the  City  University  of  New  York,  365  Fifth  Ave,  New  York,  NY  10016, 

USA; 

cPhysies  Department,  Fairfield  University,  1073  North  Benson  Road,  Fairfield,  CT  06824,  USA 


ABSTRACT 

Time  Reversal  Optical  Tomography  (TROT)  is  developed  to  locate  extended  target(s)  in  a  highly  scattering  turbid 
medium,  and  estimate  their  optical  strength  and  size.  The  approach  uses  Diffusion  Approximation  of  Radiative  Transfer 
Equation  for  light  propagation  along  with  Time  Reversal  (TR)  Multiple  Signal  Classification  (MUSIC)  scheme  for 
signal  and  noise  subspaces  for  assessment  of  target  location.  A  MUSIC  pseudo  spectrum  is  calculated  using  the 
eigenvectors  of  the  TR  matrix  T,  whose  poles  provide  target  locations.  Based  on  the  pseudo  spectrum  contours,  retrieval 
of  target  size  is  modeled  as  an  optimization  problem,  using  a  “local  contour”  method.  The  eigenvalues  of  T  are  related  to 
optical  strengths  of  targets. 

The  efficacy  of  TROT  to  obtain  location,  size,  and  optical  strength  of  one  absorptive  target,  one  scattering  target,  and 
two  absorptive  targets,  all  for  different  noise  levels  was  tested  using  simulated  data.  Target  locations  were  always 
accurately  determined.  Error  in  optical  strength  estimates  was  small  even  at  20%  noise  level.  Target  size  and  shape  were 
more  sensitive  to  noise.  Results  from  simulated  data  demonstrate  high  potential  for  application  of  TROT  in  practical 
biomedical  imaging  applications. 

Keywords:  Diffuse  optical  imaging,  time  reversal,  optical  tomography,  biomedical  imaging,  breast  cancer,  Multiple 
Signal  Classification  (MUSIC),  near-infrared  imaging 


1.  INTRODUCTION 

Time  Reversal  Optical  Tomography  (TROT)  has  been  developed  [1]  to  detect  and  locate  small  targets  embedded  in 
highly  scattering  turbid  media.  The  method  was  based  on  Diffusion  Approximation  of  the  Radiative  Transfer  Equation 
(RTE)  to  describe  light  propagation  in  a  highly  scattering  turbid  medium,  and  a  Time  Reversal  (TR)  Multiple- Signal- 
Classification  (MUSIC)  algorithm  developed  by  other  groups  in  acoustics  and  radar  applications  [2-7]. 

In  this  paper,  we  extend  the  TROT  approach  further  to  detect,  locate  and  retrieve  size  and  optical  property  information  of 
extended  targets  embedded  in  a  turbid  medium.  We  test  the  formalism  so  developed  using  simulated  data  for  a  single 
absorptive  target,  a  single  scattering  target,  and  two  absorptive  targets  assuming  different  noise  levels.  This  paper  is 
organized  as  follows.  Section  2  outlines  the  TROT  formalism  for  extended  targets.  In  section  3,  simulated  data,  TROT 
analysis  and  results  are  presented.  Section  5  serves  as  discussion  and  summary. 

2.  FORMALISM 

TROT  formalism  for  locating  small  (point-like)  targets  has  been  detailed  in  our  earlier  publication  [1],  and  may  be  used 
to  locate  extended  targets  as  well.  Here  we  present  a  brief  overview  of  the  formalism  for  completeness,  and  outline  how 
the  size  and  optical  strength  of  an  extended  target  may  be  estimated.  The  propagation  of  a  near-infrared  (NIR)  beam  of 
light  through  a  highly  scattering  turbid  medium  with  embedded  targets,  whose  optical  properties  are  different  from  that 
of  the  intervening  medium,  may  be  approximated  to  be  the  diffuse  transmission  of  light  through  background  medium  of 
uniform  optical  characteristics,  with  targets  as  perturbations.  In  the  first  order  Born  approximation,  the  perturbation  in 
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the  light  intensity  distribution  due  to  the  presence  of  the  targets  (inhomogeneities  in  optical  properties)  can  be  expressed 
using  a  data  matrix  [1]  of  the  form: 


{M  ~]  M 

lPd^^mymGs{xm,rj)  =2>„(j 

m= 1  J  m- 1 


Lm  )  ^m&s  \^m  )  ’ 


(la) 


for  absorptive  targets,  and 

M 

z  sagMmymdagsT (xm),  (ib) 

m=\a=[x,y,z } 

for  scattering  targets,  where  g5(r)  =  {G\r,  rj),  j  =  1,  . . .,  Ns}  and gd(r )  =  {Gd(rh  #•),/  =  1,  ...,Nd}  are  the  Green’s  function 
vectors  (GFVs)  associated  with  the  source  and  detector  planes,  respectively;  the  superscript  T  denotes  transpose;  G\r ,  rs) 
and  Gd(rd ,  r )  are  the  Green’s  functions  that  describe  light  propagations  in  the  background  medium  from  a  source  at  rs  to  a 
target  (inhomogeneity)  at  r  and  from  the  target  to  a  detector  at  rd,  respectively;  rm  =  8jua  (Xm  )cSVm  (rm=  SD(Xm  )cSVm  ) 
is  the  absorptive  (scattering)  optical  strength  of  the  mth  absorptive  (scattering)  target  at  with  volume  SVm;  Sjua  ( SD )  is 
the  difference  in  the  absorption  (diffusion)  coefficients  between  the  target  and  the  background  medium;  Ns,  Nd  and  M  are 
the  numbers  of  sources,  detectors  and  targets,  respectively.  It  is  assumed  the  number  of  targets  is  less  than  the  number  of 
sources  and  detectors,  M<  min(A^,  Ns);  c  is  the  light  speed  in  the  medium.  A  time  reversal  matrix  is  constructed  as  TSDDS 
=  KfK  [ TDSsd  =  {Kt)!Kt  =  K*Kr]  in  frequency  domain,  and  TSdds  =  KT K  ( TDSsd  =  KKT)  in  the  continuous  wave 
illumination.  TDSsd  and  TSdds  have  a  common  set  of  eigenvalues  {^,  j  =  1,  ...,  min(A5,  Nd)},  and  different  sets  of 
eigenvectors  {uh  i  =  1,  ...,  Nd}  and  {v7,  j  =  1,  ...,  Ns },  respectively.  The  eigenvectors  are  separated  into  signal  and  noise 
subspaces  using  an  Z-curve  method  [8]  with  an  eigenvalue  threshold  s.  For  absorptive  targets,  the  locations  are 
determined  using  the  MUSIC  pseudo  spectrum  [1] 


associated  with  the  source  plane  or  a  similar  form  for  the  detector  plane  Pd(Xp ),  or  the  product  of  these  two, 

p(xp)=p'(xp)p*(xp)’  ^ 

where  Xp  is  a  test  target  position  in  the  sample  volume.  Since  the  eigenvalues  and  eigenvectors  of  TSdds  and  TDSsd  can  be 
connected  using  singular  value  decomposition  (SVD),  i.e. 

(3) 

where  P’=  {vy}  and  &  =  {«,},  corresponding  to  £  =  diag >  c!>  are  matrices  for  the  signal  subspaces.  By 

comparing  Eq.  (3)  and  Eq.  (1),  the  target  optical  property  can  be  retrieved  by  transforming  the  eigenvalue  matrix  Z  via 

r»(Grf)-1&lU((G*)7’)'1,  (4) 

where  T  =  diag{rw,  m  =  1,  ...,  M};  Gs  =  {gs(^m)},  Gd  =  {gd(rm)}  are  matrices  including  GFVs  associated  with  the 
retrieved  target  positions.  For  scattering  targets,  the  GFVs  gd  and  gs  in  Eqs.  (2)  and  (4)  are  replaced  by  d„gd  and  dags,  a  = 
x,  y,  z,  respectively. 

An  optimal  contour  (a  surface  Q  when  plotted  in  3D)  in  the  pseudo  spectrum  in  logarithmic  scale  is  selected  to  be  the 
boundary  of  the  target(s),  via  [4] 


D.  =  arg min || K  -  Kcal  (Q)f  ,  (5) 

Q 

where  K  is  normalized  data  matrix  obtained  from  known  target  surface  in  simulation  (from  experimental  measurements 
for  real  targets)  and  Kcai(f2)  is  normalized  data  matrix  calculated  from  the  contour  of  the  pseudo  spectrum  in  logarithmic 
scale.  The  Green’s  functions  used  in  the  calculation  are  those  for  the  intervening  medium. 
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3.  SIMULATIONS  AND  ANALYSIS 


The  sample  was  taken  to  be  a  40-mm  thick  uniform  scattering  slab  with  lateral  dimension  of  80  mm  x  80  mm.  Its 
absorption  and  diffusion  coefficients  were  assumed  to  be  fia  =  0.003  mm"1  and  D  =  1/3  mm  (transport  mean  free  path,  lt  = 
1  mm),  respectively,  which  are  similar  to  the  average  value  of  those  parameters  for  human  breast  tissue.  The  index  of 
refraction  n  of  the  medium  was  taken  to  be  1.33.  The  speed  of  light  is  2.998  x  108  m/s,  or  299.8  mm/ns  in  vacuum,  and 
225.4  mm/ns  in  the  medium.  Three  simulated  datasets  were  generated  with  10-mm  diameter  spherical  targets  embedded. 
In  the  first  dataset,  an  absorptive  target  was  centered  at  (40,  40,  20)  mm.  In  the  second  dataset,  two  absorptive  targets 
were  located  at  (20,  40,  20)  mm  and  (60,  40,  20)  mm,  respectively.  In  the  third  dataset,  a  scattering  target  was  centered  at 
(40,  40,  20)  mm.  The  absorption  coefficient  of  all  the  absorptive  targets  was  set  to  be  higher  than  the  background  with 
A na  =  0.001  mm'1,  while  the  diffusion  coefficient  was  taken  to  be  the  same  as  that  of  background.  The  diffusion 
coefficient  of  the  scattering  target  was  set  to  be  lower  than  the  background  (higher  scattering  coefficient)  with  AD  =  -0.1 
mm  (lt  =  0.7  mm),  while  the  absorption  coefficient  was  taken  to  be  the  same  as  that  of  the  background.  The  volume  of  all 
targets  was  515  mm3  when  the  sample  volume  is  discretized  into  1  mm  x  1  mm  x  1  mm  voxels  in  the  forward  model. 
The  optical  strength  of  the  absorptive  targets  was  AjuacAV=  116.08  mm3/ns;  while  the  optical  strengths  of  the  scattering 
target  was  ADcAV=  -11608.1  mm5/ns.  The  incident  CW  beam  step  scanned  the  sample  at  41  x  41  grid  points  covering 
an  80  x  80  mm2  area,  with  a  step  size  of  2  mm.  Light  on  the  opposite  side  was  recorded  at  41  x  41  grid  points  covering 
the  same  area.  Additive  Gaussian  noise  of  different  noise  levels  was  added  to  the  simulated  data.  The  data  matrix  K  was 
then  obtained  using  Eq.  (1)  directly,  and  analyzed  using  TROT.  The  results  are  shown  below.  For  simplicity,  when  the 
reconstruction  result  for  one  target  is  displayed,  a  smaller  volume  of  40  mm  x  40  mm  x  40  mm  around  the  center  is 
selected.  Due  to  the  distortion  in  the  retrieved  target  shape,  target  volume  will  be  used  to  describe  the  target  size. 

3.1  One  absorptive  or  scattering  target 

The  first  20  eigenvalues  of  the  TR  matrix  for  one  absorptive  target  with  20%  added  noise  are  plotted  in  logarithmic  scale 
and  shown  in  Fig.  1(a).  The  dimension  of  the  signal  subspace  is  determined  to  be  3.  Using  Eq.  (2),  the  pseudo  spectrum 
was  calculated.  Both  of  axial  and  sagittal  views  of  the  target  using  the  pseudo  spectrum  are  plotted  in  logarithmic  scale 
in  Figs.  1(b)  and  1(c). 


Fig.  1.  (a)  Eigenvalues  plotted  in  logarithmic  scale;  (b)  and  (c)  are  the  axial  and  sagittal  views  of  the  target  using  the  pseudo 
spectrum  in  logarithmic  scale;  (d)  is  the  retrieved  target  image 

The  center  of  the  target  was  accurately  determined  to  be  (40,  40,  20)  mm.  Using  Eq.  (4),  the  optical  strength  was  found 
to  be  1 12.4  mm3/ns  when  all  3  eigenvalues  and  eigenvectors  in  the  signal  subspace  were  used,  as  shown  in  Table  1.  We 
also  found  that  if  only  the  first  eigenvalue  and  eigenvector  in  the  signal  subspace  were  used,  the  optical  strength  was 
found  to  be  112.4  mm3/ns  as  well.  From  this  observation,  we  conclude  that  the  property  information  of  the  extended 
target  is  mainly  contained  in  the  first  eigenvalue,  i.e.  the  number  of  the  eigenvalues/eigenvectors  used  is  the  same  as  the 
number  of  targets. 

Estimation  of  the  target  volume  begins  with  choosing  the  optimal  contour.  Starting  from  the  maximum  of  the  pseudo 
spectrum,  successive  contour  levels  are  selected.  When  a  contour  which  is  a  surface  Q  plotted  in  3D,  is  selected,  the 
volume  enclosed  inside  Q  is  assumed  to  be  the  target  volume,  and  then  the  data  matrix  is  calculated  and  compared  with 
the  working  (simulated  or  experimental)  data.  When  the  next  bigger  volume  is  selected,  only  the  signal  generated  by  the 
extra  volume  dV  is  calculated.  By  using  this  “local  contour”  method  [6],  the  optimal  contour  is  found  using  Eq.  (5).  The 
reconstructed  3 D  image  of  the  absorptive  target  is  displayed  in  Fig.  1(d).  The  retrieved  target  volume  is  shown  in  Table 
1  and  compared  with  the  known  volume. 

Similar  simulated  data  was  generated  and  the  subsequent  analysis  was  conducted  when  no  noise,  5%  noise,  and  10% 
noise  was  added  separately.  Similar  images  were  obtained  (not  shown).  The  target  location  and  optical  strength  of  the 
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target  was  accurately  retrieved  in  all  cases.  The  size  in  all  cases  was  also  retrieved.  The  retrieved  target  optical  strength 
and  size  are  all  shown  in  Table  1,  and  compared  to  the  known  values. 

Similar  simulations  with  one  scattering  target  were  carried  out.  The  eigenvalues  and  eigenvectors  of  the  TR  matrix,  and 
the  pseudo  spectrum  were  calculated.  The  target  location  of  the  scattering  targets  was  accurately  found  to  be  (40,  40,  20) 
mm  at  all  noise  levels  (0%,  5%,  10%  and  20%).  Axial  and  sagittal  images  of  the  scattering  target  were  obtained  using 
the  pseudo  spectrum  as  was  done  for  the  absorptive  target.  The  scattering  optical  strength  of  the  target  was  retrieved 
within  5.4%  error  and  the  size  of  targets  was  found  within  30%  error  for  all  noise  levels. 

Further  simulations  show  the  optical  strength  can  be  retrieved  even  if  the  noise  level  is  further  increased  or  the  target 
size  and  location  vary;  however,  the  error  in  the  retrieved  target  size  may  increase  dramatically  if  the  target  size  further 
increases.  The  accuracy  of  the  retrieved  target  size  and  optical  property  seems  not  to  well  correlate  with  the  noise  level. 


Table  1.  Optical  strength  and  size  of  an  absorptive  target  and  a  scattering  target  (10-mm  diameter). 


Noise  Level 
(%) 

Absorptive  Target 

Scattering  Target 

Optical  Strength 

Size 

Optical  Strength 

Size 

Retrieved 

(mm3/ns) 

Error 

(%) 

Retrieved 

(mm3) 

Error 

(%) 

Retrieved 

(mm5/ns) 

Error 

(%) 

Retrieved 

(mm3) 

Error 

(%) 

0 

112.4 

3.2 

478 

7.0 

-10985.4 

5.4 

413 

19.8 

5 

112.4 

3.2 

484 

6.0 

-10985.8 

5.4 

535 

3.9 

10 

112.4 

3.2 

370 

28.2 

-10985.2 

5.4 

665 

29.1 

20 

112.4 

3.2 

385 

25.2 

-10987.2 

5.3 

576 

11.8 

*  Known  values:  volume:  515  mm3,  absorptive  optical  strength:  1 16.08  mm3/ns,  scattering  optical  strength:  -1 1608.1  mm5/ns 


3.2  Two  absorptive  targets 


Next  we  considered  the  case  of  two  spherical  absorptive  targets  (diameter  10  mm)  embedded  in  the  medium,  with  a 
center-to-center  separation  of  40  mm.  Different  added  noise  levels:  0%,  5%,  10%,  20%,  and  100%  were  considered.  The 
eigenvalues  and  eigenvectors  of  the  TR  matrix,  and  the  pseudo  spectra  were  calculated.  The  target  locations  of  two 
targets  were  accurately  found  to  be  (20,  40,  20)  mm  and  (60,  40,  20)  mm  at  all  noise  levels.  Typical  axial  and  sagittal 
images  of  the  targets  generated  using  the  pseudo  spectra  for  a  noise  level  of  20%  are  displayed  in  Fig.  2.  Similar  images 
were  obtained  for  other  noise  levels.  The  optical  strength  and  size  of  both  targets  were  found.  The  optical  strength  of 
each  target  was  calculated  within  3.3%  error  at  all  noise  levels.  The  size  of  targets  was  found  within  13.4%  error  when 
no  noise  was  present,  and  varied  with  noise  added.  The  retrieved  optical  strength  and  size  of  the  targets  are  listed  in 
Table  2  for  different  noise  levels.  Since  both  targets  have  the  same  optical  property  and  size  in  the  simulated  data  and  the 
retrieved  data,  only  one  set  of  values  is  shown  in  the  table. 


sagittal  (x  =  20  mm)  sagittal  (x  =  60  mm) 


20  40  60  80  20  40  60  80 

y  (mm)  y  (mm) 


(c) 


(d) 


Fig.  2.  (a)  Eigenvalues  plotted  in  logarithmic  scale;  (b),  (c)  and  (d)  are  the  axial  and  sagittal  views  of  the  targets  using 
pseudo  spectrum  in  logarithmic  scale;  (e)  is  the  retrieved  target  image.  The  added  noise  level  was  20%. 

Further  simulations  (not  shown  here)  indicate  that  retrieved  target  size  was  more  accurate  when  small  targets  were 
involved.  Simulations  with  one  target  at  different  locations  and  two  targets  with  different  separations  were  also  tested  for 
absorptive  and  scattering  targets.  The  target  locations  were  accurately  retrieved  in  all  of  these  cases.  The  retrieved 
optical  strength  in  all  cases  was  much  less  sensitive  to  the  target  location  and  size  than  the  retrieved  target  size.  As 
expected,  it  was  more  challenging  to  retrieve  the  optical  strength  and  size  of  scattering  targets  than  absorptive  targets. 

Further  simulations  also  showed  closer  distance  between  two  targets  made  target  size  retrieval  more  difficult,  because  of 
the  cross  talk  between  the  two  targets  showing  up  in  the  contour  of  the  pseudo  spectrum  (logarithmic  scale). 
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Table  2.  Optical  strength  and  size  of  two  absorptive  targets  (10-mm  diameter). 


Noise  Level 
(%) 

Optical  Strength 

Size  | 

Retrieved 

(mm3/ns) 

Error 

(%) 

Retrieved 

(mm3) 

Error 

(%) 

0 

112.3 

3.3 

584 

13.4 

5 

112.3 

3.3 

295 

42.7 

10 

112.3 

3.3 

368 

28.5 

20 

112.3 

3.3 

487 

5.4 

100 

112.2 

3.3 

663 

28.7 

*  Know  values:  volume:  515  mm3,  optical  strength:  116.08  mmVns 


4.  SUMMARY  AND  DISCUSSION 

Time  reversal  optical  tomography  (TROT)  was  further  developed  to  deal  with  extended  targets.  The  center  position  of 
the  target(s)  is  determined  rather  accurately  for  both  absorptive  and  scattering  targets.  It  is  found  that  the  optical  strength 
(absorption  or  scattering)  can  be  retrieved  for  different  target  size,  target  location  and  noise  level.  However,  it  is  much 
more  challenging  to  retrieve  the  target  size.  The  retrieved  target  size  is  determined  by  the  details  of  the  pseudo  spectrum. 
There  seems  to  be  a  lack  of  well  defined  correlation  between  the  noise  level  and  the  size  and  optical  strength  of  the 
targets,  which  needs  to  be  understood.  The  efficacy  of  the  approach  will  be  further  tested  using  experimental  data. 
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Appendix  3 


OPTICAL  COHERENCE  TOMOGRAPHY 

Time  Reversal 

Optical  Tomography 

There  has  been  a  surge  of  interest  in 
diffuse  optical  tomography  (DOT) 
that  uses  near-infrared  light  to  detect, 
localize  and  diagnose  maladies  such  as 
breast  cancer  and  brain  injury.1  Scattering 
and  light  attenuation  limit  the  resolution 
and  accuracy  of  DOT  methods  that  use 
small  differences  in  optical  properties  to 
distinguish  lesions  from  normal  tissue. 
Researchers  need  a  DOT  approach  that  can, 
for  example,  quickly  reconstruct  images 
to  detect  and  map  tumors  at  early  growth 
stages  and  determine  if  they  are  malignant 
or  benign.2 

Time  reversal  optical  tomography 
(TROT)  extends  the  use  of  TR  imaging  and 
a  subspace-based  method  of  multiple  signal 
classification  (MUSIC)  from  acoustic  and 
radar  imaging  to  optical  imaging.3,4’5  TROT 
uses  a  multisource  illumination  and  multi¬ 
detector  signal  acquisition  scheme  to  acquire 
multiple  angular  views  of  the  sample. 

The  perturbation  in  light  intensity 
distribution  due  to  the  targets  is  extracted 
from  the  data  and  organized  in  a  matrix  K. 

The  leading  eigenvalues  of  the  TR  matrix, 

T  =  correspond  to  the  targets  whose 
locations  are  determined  using  MUSIC, 
along  with  Green’s  functions  for  light 
propagation  in  the  sample. 

We  first  tested  the  efficacy  of  TROT  using 
a  60-mm-thick  slab  of  Intralipid-20  percent 
suspension  in  water  and  9-mm  diameter 
glass  spheres  as  absorptive  or  scattering 
targets.  We  filled  the  glass  spheres  with 
ink  dissolved  in  the  suspension  to  provide 
absorptive  targets,  and  with  a  higher  con¬ 
centration  of  Intralipid  to  provide  scattering 
targets.  We  chose  the  optical  properties  and 
size  of  the  sample  and  targets  to  emulate 
average  values  for  breast  tissue  and  small 


(a)  F  =  signal  transmitting  narrow-band  filter;  TS  =  translation  stage;  CCD  =  charge  cou¬ 
pled  device;  and  PC  =  computer.  Continuous  wave  790-nm  diode  laser  light  illuminates  the 
front  of  the  sample  cell.  Diffusely  transmitted  light  from  the  opposite  face  is  collected  by 

a  camera  lens  through  F  and  sensed  by  a  CCD  camera.  The  sample  cell  is  step-scanned 
across  the  laser  beam  in  a  2-D  x-y  array  of  grid  points  using  the  computer-controlled  TS. 

(b)  A  TROT-generated  pseudo  image  of  two  absorptive  targets  at  z  =  30.5  mm  plane  when 
the  targets  are  separated  by  27.6  mm. 


tumors.  We  found  that  TROT  could 


Researchers 


retrieve  the  location  of  a  single  target 
with  millimeter  accuracy  and  resolve  two 
targets  when  their  adjacent  surfaces  were 
only  4-mm  apart. 

Another  experiment  involved  a 
realistic  breast  model  composed  of  ex 
vivo  breast  tissue  with  two  pieces  of 
embedded  tumors;  TROT  accurately 
located  the  positions  of  both  the  tumors. 
We  have  extended  TROT  for  locating 
fluorescent  targets. 

TROT  is  non-iterative  and  faster  than 
other  iterative  DOT  approaches.  It  is 
particularly  suited  for  detecting  point¬ 
like  targets.  MSEI 
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